- Email: [email protected]
Least Square Solution of Non Linear Problems in Fluid Dynamics
G.M. de La Penha. L.A. Medeiros (eds.) Contemporary Developments i n Continuum Mechanics and P a r t i a l D i f f e r e n t i a l Equations @North-Holland P u b l i s h i n g Company (1978)
LEA.I;T SQUARE SOLUTION O F NON L I N E A R PROBLEMS
I N F L U I D DYNAMICS
R.
GLOWINSKI* and
Universitb de P a r i s Iria-Laboria,
0.
PIRONNEAU
6 and
Iria-Laboria, 78150 le Chesnay
France
France
1. I n t r o d u c t i o n . Many problems i n Mathematical P h y s i c s a r e of
the
following forms, e i t h e r A(u)u = f,
+ B(u)
Au
In
(l.l), A ( v )
i s in (1.2)
v
i s for
=
f.
g i v e n a l_i _ n e_a_ r . o p e r a t o r and
B
a non l i n___ e a r operator.
This suggests obviously the following algorithms: F o r s o l v i n g (1.1)
(1.3)
uo
then f o r
n
2
0
we o b t a i n
u n+l
given from
u
n
by
I
A(un)un+' = f ,
For s o l v i n g ( 1 . 2 )
we may use an a l g o r i t h m s i m i l a r t o
(1.5)
(1.4)1 i s r e p l a c e d by
*
except t h a t
(1.3)-
~.
V i s i t i n g P r o f e s s o r IJF'RJ,
R i o d e J a n e i r o , July-August
1977.
172
R. GLOWINSKI and 0. P I R O N N E A U
There are some applications in which (l.3), (resp. ( 1 . 3 ) , ( 1 . 4 ) 2 , (resp. (1.2)).
-
(1.5))
( 1 . 4 ) 1’ (1.5)
converges to a solution o f (1.1)
Unfortunately these are important applications
mathematically modelled by equations like (1.1) or (1.2)
-
f o r which the above algorithms can not be used, therefore one
has to use more sophisticated methods.
In some cases (cf., e.g., Glowinski-Marocco Ell) these problems may be reduced to a saddle-point problem by introduc-. ~
ing an artificia-1 constraint and the corresponding Lagran* multiplier. - __ -..~
In some other situation such an approach is not suitable (by lack o f monotonicity usually) and it may be convenient to use aleast square method.
A typical least square formulation o f problem (1.1) is Min J(v)
(1.6)
V
where
(1.7) In
(1.7),
]I*]]
denotes a suitable norm (for more details
about this last point, see Remark 1.2 below) and
J(v)
the
solution o f the linear problem
We obtain similarly a least square formulation of (1.2) by replacing (1.8)1 by
(1*8)2
Ay = f-B(v).
LEAST SQUARE S O L U T I O N O F NON L I N E A R PROBLEMS
Remark _ _ 1~ .1:
I t a p p e a r s c l e a r l y from t h e above r e l a t i o n s t h a t
(1.6), (l.7), s t r u c t u r e of
(l.7),
( 1 . 8 ) l and ( 1 . 6 ) ,
( l . 8 ) 2 have t h e
p t i m a l c o n t r o l problem i n which t h e c ontrol an o ._ ~~
variable
(i.e.
t h e independent v a r i a b l e i n t h e problem) i s
the s t a t e variable i s
y
and t h e c o-s t-f-u n-c.t_ i -o. n is -
v,
J.
T h i s p a r t i c u l a r s t r u c t u r e l e a d s v e r y n a t u r a l l y t o compute t h e derivative o f
111 , [ 2 ] )
J
v i a an a d j o i n t s t a t e e q u a t i o n ( s e e Lions
f o r t h e o p t i m a l c o n t r o l o f s y s t e m s m o d e l l e d by
p a r t i a l d i f f e r e n t i a l equations). T h i s k i n d o f methodology h a s b e e n a d v o c a t e d by Cea-Geymonat
[l] f o r s o l v i n g non monotone, non l i n e a r t w o - p o i n t b o u n d a r y v a l u e problems l i k e
1
(1.9)
+
-u”
cp(u) = f
over
10,1[
u ( 0 ) = u ( 1 ) = 0.
I n B e g i s - G l o w i n s k i [1] a s i m i l a r a p p r o a c h h a s b e e n u s e d f o r s o l v i n g f r e e boundary problems
l i k e those r e l o t e d t o flows o f
l i q u i d s i n porous media. I n a n u n p u b l i s h e d work Derviauc-Glowinski h a v e a l s o u s e d t h i s t y p e o f methods t o compute t h e maximal s o l u t i o n of
I
(1.10)
where
>
-u
= Xe
U
over
10,1[
”
u(0) = u ( 1 ) = 0
0.
In t h e above r e f e r e n c e s t h e o p t i m i z a t i o n i s a c h i e v e d u s i n g e s s e n t i a l l y -~ descent o r steepest descent algorithms. ~
~~
~~~
More
~~~~
r e c e n t l y G l o w i n s k i and P i r o n n e a u o b s e r v e d t h a t t h e u s e of a conjugate gradient algorithm,
t o achieve the minimization,
c a n d r a s t i c a l l y improve t h e s p e e d o f c o n v e r g e n c e ,
specially
17 4
R. GLOWINSKI and 0. PIRONNEAU
in the solution of some hard problems like those to be considered in the following sections.
J
Remark 1.2: The choice of a convenient cost function
is a
very important matter if one wishes to use the above rnethodology.
Let us mention some possibilities we have found very
attractive since they simplify the adjoint state equation and lead to good convergence properties.
In the case of the problem (1.1) assume that for A(v)
E e(V,V’)
dual space. if
A(v)
where
(
where
is an Hilbert space and
V
Assume also thak
A(v)
v
given, V’
is self-adjoint. .
is _ strongly - a natural choice for _ _ _ V-elliptic
,*)
denotes the cu-ality -~_ pairing ~ _ . _o _ f _V’
Similarly i f in (1.2),
A
and
its Then
J
is
V.
is strongly V-elliptic it is
convenient to use
In the following sections we shall see how the above principles can be applied to the numerical solution of two important problems i n Fluid Mechanics, namely: *
The transonic flow of an isentropic, compressible perfect ~
fluid. ’ Navier-Stokes equations for incompressible fluids.
Numerical results will be presented in order to illustrate the possibilities of the above methods.
To conclude this report some indication will be given on the minimization of non quadratic functionals by conjugate gradient algorithms.
LEAST S Q U A R E S O L U T I O N OF NON L I N E A R PRODLEMS
175
2. Numerical solution-of transonic flow problems. . __ ~
2.1
Generalities on transonic flows.
The theoretical and numerical studies of transonic flows for -_ perfect fluids have always been very important,
But these
problems have become even more important these last years in relation with the design and development of large subsonic aircrafts. These transonic problem are very difficult, theoretically and mumerically, for the following reasons: They are n o n linear. Shocks may exist in the flow. An Entropy Condition is needed, in a way o r another, to avoid non physical solutions. They are mixed problems i n the sense that they are el_l.ipt-icin the Eubsgsic_ part of the flow, LyEerbolic i n the supersonic part of the flow. From a theoretical point of view let us mention the work of Morawetz [l].
At the present moment the numerical methods
which are the more commonly used have originated from Murman-Cole [l]
and we shall mention B a u e r - G a r a b e d i a n - K o r n r l ] ,
Bauer-Garabedian-Korn,Jameson
11
,
Jameson [ 11 ,[ 21 ,[ 31 ,[ 41
and the bibliographies there i n (see also Hewitt-Illingworth and Co-editors [ 13 ) . These above numerical methods use the key idea of Murman and Cole which consists to use a finite difference scheme, centered i n the subsonic-part of the flow, kackward (in the direction of the flow) in the supersonic part.
The switching
1 76
R.
GLOWINSRI a n d 0 .
between t h e s e two-schemes
PIRONNEAU
i s a u t o m a t i c a l l y done v i a a
t r u n c a t i o n o p e r a t o r only a c t i v e i n the supersanic p a r t of flow ( s e e Jameson,
loc.
cit.,
f o r more d e t a i l s ) .
We s h a l l d e s c r i b e h e r e a d i f f e r e n t a p p r o a c h ,
which seems v e r y
c o n v e n i e n t f o r comput i n g f l o w s p a s t p r o f i l e s ,
or f l o w s i n n o z z l e s .
infinity,
the
subsonic a t
I n o u r method t h e t r a n s o n i c
f l o w problem i s f o r m u l a t e d a s a non l i n e a r l e a s t s q u a r e problem,
control w h i c h i s i n t u r n i n t e r p r e t e d a s a n o p t i m a l~-
p r o b l -e m .
Then f i n i t e e l e m e n t s a r e u s e d t o a p p r o x i m a t e t h e Since t h e e n t r o p y~-c o n d i t i o n i s f o r m u l a t e d by .____
above problems.
~
a l i_ n e_ ar _ ine q u a l i t y c~-o n s t . ra i n t , a c o n v e n i e n t method t o h a n d l e .
i t i s t o u s e p e n a l t y a n d / o r d u a l i t y methods ( u s i n g a n augmented l a g r a n g i a n i f p e n a l t y and d u a l i t y a r e co mb in ed ) . Duality techniques a r e not
c o n s i d e r e d i n t h i s r e p o r t and w e
r e f e r t o Glowinski-Pironneau [ d e s c r i p t i o n of
2.2 2.2.1
11,
Glowinski
[13 for a
them.
M a t h e m a t i c a l m o d e l s f o r t r a n s__ o n __ i c f.l_ o ws. . .. B..-a s i c a s s u m p t i o n s .
We a s s u m e t h a t t h e f l u i d u n d e r c o n s i d e r a t i o n i s p e r f e c t a n d compressible,
and t h a t t h e flow of
a n d i r r o t a t i o-n a l ( i . e . ,
potential).
-
such a f l u i d is isentropic These assumptions a r e n o t
t r u e i n g e n e r a l s i n c e t h r o u g h a s h o c k wave t h e r e i s a v a r i a t i o n of e n t r o p y , and a n i r r o t a t i o n a l f l o w becomes rotational
(i.e.
t h e v a l i d i t y of i n t h e c a s e of
t h e r e i s c r e a t i o n o f v o_ r t_ i c i_ ty).
Therefore
t h e model t o f o l l o w i s assumed t o b e c o r r e c t "weak s h o cks". -.
177
LEAST SQUARE SOLUTION O F N O N LINEAR PROBLEMS
I n t h e c a s e o f a f l o w s p a s t a s- h a_ r p~ p_r o_f ._ i l e w e s h a l l suppose t h a t t h e r e i s no wake b e h i n d
2.2.2
E--q u a -~ t i o n s of
the tr a- _ i l ~i n g edge.
t h e flow.
~
n
Let
be t h e domain o f
r
t h e f l o w and
i t s boundary;
then
we s h a l l assume t h a t t h e f l o w i s m o d e l l e d by
-
(2.1)
v-((1 -
,
1
, 2
y-1
where i n ( 2 . 1 ) .
- 0 -
C,
-
y
the flow v e l o c i t y ,
i s t h e c.r i t i c a l v e l o c i t y , ~~
i s t h e r. __ a t i o of
W e have *
09
i s the flow p otential, _____
(y = 1 - 4 f o r
specific heats _ _ I
4.
t o add t o ( 2 . 1 )
Boundary c o n d i t i o n s ( o f D i r i c h l e t a n d / o r
Neumann t y p e , f o r
instance), *
-
A K u t t a - J o u k o v-k y
condition i n the case of
a l i f t i n g body ( s e e L a n d a u - L i f c h i t z
the flow past
[ l , Sec.
46]),
An e n t r o p y c o n d i t i o_ n . i n o r d e r t o e l i m i n a t e t h e non p h y s i c a l .s o l u t i o n s of
(1.1); t h i s p o i n t
w i l l be d i s c u s s e d i n S e c .
2.2.3.
:
Remark 2 . 1 :
I t c a n h a p p e n t h a t on some p a r t of
0
have
and
ness; -
t o be g i v e n s _ i m_ u l t a n e o u s- __ ly t o i n s u r e unique-
i t is, f o r i n s t a n c e ,
of F i g .
2.1 i f
~~
t h e c a s e of
'4
g i v e n on
t h e divergent nozzle
the velocity a t the entrance i s supersonic.
T y p i c a l boundary c o n d i t i o n s a r e 5-n
t h e boundary,
r4, rl, r 2 .
If
@
g i v e n on
r l , Fg
and
the velocity a t the entrance
178 (i.e.
R.
I-,)
GLOWINSKI and 0. P I R O N N E A U
is subsonic we require ___ less boundary conditions.
Figure 2.1
Figure 2.2 Remark 2.2:
In the case o f t h e f l o w past a multibody profile
(like i n Fig. 2.2)
each
body requires a Kutta-Jourkowsky
condition.
2.2.3
Formulation of the entropy condition.
It follows f r o m Landau-Lifchitz [ l , ch.91 condition c a n be formulated as follows:
that the entropy
1 79
LEAST SQUARE SOLUTION OF NON LINEAR PROBLEMS
1 In the direction
I
-!
(2.2)
o f the flow, One- cannot have a - ___ ~~
subsonic-supersonic transition through
a_;t?ck.
F o r a one-dimensional . . flow, (2.2) implies ~
*<
2 d
(2.3)
dx
3-
2
+m,
is a measure bounded from_ above; _ ~ _ _weak (and more dx rigorous) formulation of (2.3) are: i.e.
-~
There exists a constant ~~
~
M,
such that-either -
where
(We recall that fJ
(R) =
{cp E Cm ( f i )
,
cp
has a compact support in
n] . )
In the case o f a t ~ oo r three dimensional flow we shall suppose that (2.2) can be formulated by
(2.9)
190
R.
GLOWINSKI a n d 0 . P I R O N N E A U
T h e n u m e r i c a l r e s u l t s o b t a i n e d f o r t__ wo and _ t r e_e -_ d i_ m e~n s_ i o_ n a_ l flows u s i n g d i s c r e t e analogous of above formulation Remark 2 2 :
(2.7) seanto justify the
of t h e entropy condition.
The c o n d i t i o n ( 2 . 7 )
i s a c t u a l l y a p a r t i c u l a r case
O f
(w)+ E LPP)
(2.10) where,
assuming t h a t
is a m a s_ u_ re, _e_
( ~ c p ) + = positive p a r t of
(2.11) Then
Acp
(2.7)
corresponds t o
f r o m a computational point
p =
+m;
~cp.
another useful choice i s ,
of view,
p=2
i.c.
(2.12)
To c o n c l u d e t h i s S e c .
2.2
a t t h e moment,
l e t u s mention t h a t ,
t h e t r a n s o n i c f l o w problem i s s t i l l open from a t h e o r e t i c a l ~
~~
of view, even i n t h e simpler c a s e of t h e so-called
point
s m a l l d i s t u r b a n c_e _ c c l u a t i o n ~~
details).
( s e e J a m e s o n [ 11 ,[
41
f o r more
T h i s l a c k o f e x i s t e n c e and u n i q u e n e s s r e s u l t s h a s
v e r y l i k e l y i t s c o u n t e r p a r t i n t h e n u m e r i c a l methodology.
2.3
L - east
s q u a r e f o r m u l a t i on. ~
~~
R edu5;ion _ _t o~ a n o p t i m a l -~ ~~
~
~
~~
c o n t r o l problem. If
w e s u p p o s e t h a t t h e d e n s i t y o f t h e f l u i d i s o n e when t h e
velocity is zero,
then the coefficient
a s the density of
the
(2.13)
fluid.
of
V@
i n (2.1)
We s h a l l u s e t h e n o t a t i o n 1
lvcpl
2 Y - l
)
.
appm
181
LEAST SQUARE SOLUTION OF NON L I N E A R PROBLEXS
The i d e a o f t h e method t o f o l l o h i s t o d ecoup_ le t h e _ density ~ and t h e p o t e n t i a l those described
4
potential
0
and
-
via a _ l e_a s -~ t s_ q u_ a r_ e~ f o_ rm u_ l a.t_____ i o_ n like _ _
@
To d o s o we i n t r o d u c e a new
1.
i n Sec.
-
the c -o .n t r o l potential ~
by m i n i m i z i n g a c o n v e n i e n t
least
square c o s t function.
W e may u s e f o r i n s t a n c e t h e f o l l o w i n g f o r m u l a t i o n formulations s e e Glowinski-Pironneau pironneau
5
and t r y t o r e c o u p l e
(for other
[ 2 1 , Gowinski-Periaux-
[ 11 )
(2.14) where
(2.15) i
In
(2.14
,
a
'' c o n v e n i e n t l y ''
p ep(x)) = 0
y =
that for a i r , appears t h a t
via the state
p l u s boundary c o n d i t i o n s f o r
t h e parameter
a convex s e t Since
5
i s a f u n c t i o n of
in
i n the
is either cho s e n
.
implies
1
y+1 2
(--)
Y-1
t r a n s o n i c range
=
1
= (Fl) --
,,/6 = 2 . 4 5 , =
(say
X
and
1 Y+1 2
I Vcp(x) 1
i f and o n l y i f
1.4
or
0
I-.
on
d,
is
and
C,
it
lVCpl
1.5
c*)
w e have
(2.16)
o
< 6 s p(cp(x)) s 1
I t f o l l o w s from ( 2 . 1 6 )
that
a.e.
on
0.
( 2 . 1 5 ) i s an e l l i p t i c p r o b l e m f o r
a p p r o p r i a t e boundary c o n d i t i o n s . l i f t i n g b o d i e s , Kutta-Joukowsky
I n t h e case of flows p a s t conditions are a l s o required
physical t o h a v e , w i t h t h e o t h e r b o u n d a r y c o n d i t i o n s , a -__ s o l u t i o n o f problem ( 2 . 1 5 )
(modulo a c o n s t a n t i f
Neumann c o n d i t i o n s on t h e b o u n d a r y ) .
one o n l y h a s
R.
182
2.4:
Remark
If
GLOWINSKI and 0 . P I R O N N E A U
i n t h e o r i g i n a l problem,
0
r,
be s i m u l t a n e o u s l y p r e s c r i b e d on some p a r t of approach,
5
w i t h t h e two p o t e n t i a l
and
'@ have t o an
and
a,
t h e above
i s very
convenient s i n c e t h e boundary c o n d i t i o n s can be s p l i t t e d b e t ween
@
4.
and
However i f
boundary c o n d i t i o n s f o r t a k e a c c o u n t of
0
one w i s h e s t o u s e t h e same and
4
i t i s always p o s s i b l e t o
t h e e x t r a boundary c o n d i t i o n s
(assumed t o be
(2.14) a
o f D i r i c h l e t t y p e ) by a d d i n g t o t h e c o s t f u n c t i o n I
quantity proportional t o e i t h e r
( I
' rd
15-@,1
(9-0,1
2
br,
or
'd (or t o a l i n e a r c o m b i n a t i o n o f b o t h ) , where
2
i s the part
of
on which one r e q u i r e s
A s i m i l a r i d e a i s used i n Begis-Glowinski
@Ird
=
rd
Od.
[l] t o s o l v e some
f r e e boundary p r o b l e m s . Remark 2.5: variant
To s t a t e t h e e n t r o p y c__ o n d-i t_i_ o n ( 2 . 7 )
4
( 2 . 1 2 ) ) we h a v e t h e c h o i c e between
(or its
and
F
(actud-
l y we can a l s o u s e t h e s e two p o t e n t i a l s ~ i m y l t a n e o u ~ l y ) I. f one u s e s
E
( r e s p . @ ) we have a c o n t r o l c o-n s t r a_i_n t
(resp.
a s t a t e c-onstr-int). Remark 2 . 6 :
W e observe t h a t for t h e c l a s s o f flows under
c o n s i d e r a t i o n we h a v e
T h i s o b s e r v a t i o n s u g g e s t s t h a t t h e convex
X
( 2 . 1 4 ) h a s t o be t a k e n a s a convex s u b s e t of
occuring i n
Wl'"(S2).
We
o b s e r v e a l s o t h a t t o s t a y i n t h e t r a n s o n i c r a n g e i t may b e convenient t o intr oduce t h e following c o n s t r a i n t s
or
183
LEAST SQUARE S O L U T I O N OF NON LINEAR PROBLEMS
(2.18) Actually the computations we have done prove that for a physically well-posed transonic problem it is not necessary ~~
to introduce
Remark
~
(2.17), (2.18).
2.7: If the transonic problem has a solution and if
X
is "large enough" then the least square problem has a solution such that the cost function is equal to zero.
This last
property gives us means to check the quality of the computed solution.
2.4
Finite element approximations. __
We assume that
0c
_. Three-dimensional calculations will
2
!R
be reparted i n a forthcoming paper.
2.4.1
Generalities.
The above control problem will be approximated by finite elements.
Compared to finite differences, finite elements
give us the possibility o f handling problems posed on rather complicated geometrical domains.
Moreover the variational . .~
framework of finite element formulations is very appropriate to the problem under consideration and to the methodology we use to solve it.
It will be in particular fairly easy to
approximate the weak formulations of the entropy condition (2.7) (and also the alternate entropy condition (2.12)).
If
R
is unbounded, it is replaced by a bounded domain
still denoted by
-
as large as possible.
-
To approximate
the above continuous problem, we introduce a standard tri-
184
R. GLOWINSKI and 0 . PIRONNEAU
Zh
angulat ion elements).
of
61
(we can also use quadrilateral finite -
7
Then we approximate the functions
and
$I
by
piecewise polynomial functions belonging to the following subspace of
with
H1(n)
(and
Wl’”(n))
Pk = the space o f polynomials of degree
Sk.
Remark 2.8: In the case of a lifting body, to take into account the Kutta-Joukowsky condition, one usually introduces -~ .. (see Fig. 2.3) an arc
y
between the trailing edge o f the
body and the external boundary.
Thus arc
constant jump, a priori unknown of
a@
requires the continuity of
Y.
Since
@
-
an
@
(and
is discontinuous along
used any more, but introducing
.
.
0
and
y
supports a
<,
ae z) at y,
but one the crossing o f
Hl(61)
c a n not be
\
n = n - y one can use
A.
HL(n)
and define
Vh
over a triangulation of
Results of computations for such profiles are given in
Sec. 2.6;
however the method f o r t r e a t i n g t h e Kutta-Joukowsky
condition will not be treated in this report since they are not specific of transonic flows. In order to simplify the presentation o f our method we shall assume i n the sequel that one works directly over
61.
LEAST SQUARF S O L U T I O \
@rU O U
Figure
Remark 2 . 9 :
If
@
and
?
2.3
have t o s a t i s f y o n l y N e u m-a r i r i
b-o u n d a r y_ _c o n d i t -~ i o n s , t h e s u b s p a c e of ~
Vh
itself.
If
on t h e c o n t r a r )
i$
Vh
t o h e u s e d w i l l be
and/or
haxe to satisfi
F
D i r i c h l e t boundary c o n d i t i o n s over s o m e p a r t of s h a l l have t o use subspaccs o f Remark 2 . 1 0 : domain o f
R
Vh,
,
",
t h e n we
s t r i c t l y included _-___
We have t a c i t l y assumed t h a t 2
185
LIhFAR l'ROI3LEVS
R
in V
h'
i s a polygonal
or has b e e n a p p r o x i m a t e d b y s u c h a d o m a i n .
Hoizever i n t h e c a s e o f a curx-ed b o u n d a r y i t i s a l w a y s r > o s s i b l e t o use
( a t s o m e e x t r a computational c o s t ) curved f i n i t e
elements
t-11).
(see,
for i n s t a n c e , S t r a r i g - F i x [I], C i a r l e t - R a v i a r t
R. GLOWINSKI and 0. PIRONNEAU
185
2.4.2
o f the state equation . and . . o f_ _the cost Approximation __. ~~
function. T o simplify the presention, we shall assume that we only have
Neumann .boundary cond>t-ions, i.e.
aa@ n
(2.20)
= g
on
7
(it is the case for the important application of flows subsonic at infinity, past profiles). F
sections that
and
0
It follows from the above
will satisfy the same boundary ~
conditions.
We shall also assume that if
then the corresponding value of (2.21)
p g
p
dr =
g(x)
f 0, x E
r,
is known, and that 0.
'iThen the state equation (2.15) h a s the following variational
(2.22)
which is approximated by
r
i
4
(2.23) where @
and
P(ch)v@h'v"h
f 'I-
(pg)h q h
[
Oh vh' ( ~ g )is~ a Convenient approximation of @h
dr
pg.
qhevh,
Since
are only defined modulo a constant we shall
prescribe the values of some point of ed by
dx =
r.
@
and
cph
(and
5
and
5,)
at
The cost function in (2.14) is approximat-
187
LEAST SQUARE S Q L U T I O N O F PJON L I N E A R PROULEFIS
If one
R e m a r k 2.11:
(k=l),
uses
the piecewise linear approximation
then the integrals occuring i n ( 2 . 2 3 ) ,
to compute since property for
05,
p(ch)
(2.24)
piecewise constant implies the same (from (2.1'3)).
If
k=2
a ~ numerical _
integration procedure has to be used in ( 2 . 2 3 ) ,
a
are easy
_
and also in
= 1.
(2.24)
if
2.4.3
Approximation of .the entropy condition. -~ _. ~~
~
~~~
~
T o avoid non physical shocks (i.e. shocks for which the
thermodynamical entropy condition is n o t satisfied) we use the formulations o f Sec. 2 . 2 . 3 .
We still assume that w e o n l y
have Neumanii b o u n d a r y conditions.
We use the notation of the continuous oroblem. the condition ( 2 . 1 2 ) (2.14)
If
one uses
we may "regularize" the cost function
by adding to it, with
C
> 0, either
(2.25)
or (2.26)
E
[ '*
l ( h E ; ) + 1 2 dx
( o r a linear combination of both).
In
(2.25),
(2.26),
is
a "small" parameter.
We can make this approach more sophisticated by using instead of ( 2 . 2 5 ) ,
(2.26),
regularization functionals like
188
R.
GLOWINSKI and 0 . PIRONNEAU
0.
p o s s i b l y e q u a l t o z e r o o v e r some p a r t of
To u s e t h e above methodology f o r t h e a p p r o x i m a t e problem i t i s n e c e s s a r y t o have a n a p p r o x i m a t i o n o f W e may u s e a n a p p r o x i m a t i o n a p p r o x i m a t i o n of
C i a r l e t - R a v i a r t [ 21
,
[
(2.28)
73 ) :
i s s u f f i c i e n t l y smooth, t h e n from
we h a v e
Ai/lcpdx =
'n
(2
i
epdP
-
(2.29)
I
{R
an a p p r o x i m a t i o n of
E H1(n).
-
i
'R
Ah
$I,
0Bh'VCPh dx
A@
of
cph E Vh,
bL
I n (2.29),
t o define
the function
of
r,
g
of
gh
Remark 2.13: To o b t a i n
(2.20).
t h i s can be d o n e . from ( 2 . 2 9 ) we have t o s o l v e a
Ah@,
l-~ i n e a r s ys t e m t h e m a t r i x o f which i s _s_yrnmetric, p o s i t i v e ~~
sparse, but p o t d i a g o n a l
approximation of t h e o p e r a t o r approximate
is
t h e above i d e a i s s l i g h t l y more
complicated t o a p p l y , b u t
definite,
as
we a l s o h a v e D i r i c h l e t boundary c o n d i t i o n s
If
o v e r some p a r t
cp
E Vh*
We u s e t h e same method
Remark 2 . 1 2 :
ghCPhc
=
'h'hcphdx
AhOh
+
V w - v c p dx
F
U s i n g t h i s i d e a we d e f i n e a n a p p r o x i m a t i o n €allows
[21,
e q u a t i o n ( s e e Glowinski
-
Glowinski-Pironneau [
$
L e t u s assume t h a t Grgeq-for5qJ-a-
s u g g e s t e d by F i x e d -f i n i t e e_l e_m_e n t
t h e b i h a r m o n~ic
~
AS).
(resp.
I
n
4 0 cp
h h h
dx
I).
( t h i s m a t r i x i s an If
k=l
one c a n
u s i n g t h e two-dimensional,
189
LEAST SQUARE SOLUTION O F N O N L I N E A R I'RODLEMS
t____ r a p e-z o i d a l n u m e r i c a l i n t e g r a t i o - n m e t h o d .
D o i n g s o we o b t a i n
~
s o l v i n g a l i n e a r s y s t e m w i t h a _d_i a g o n a l m a t r i x .
by
ahQh
If
t h e a b o v e r e g u l a r i z a t i o n method i s t e c h n i c a l l y more
k=2
complicated t o use.
A
Once
h a s b e e n a p p r o x i m a t e d we add t o t h e c o s t f u n c t i o n
(see ( 2 . 2 4 ) ) the functional
2
I (Ahah)+]
('. 3')
dx
(rcsp'
I
Eh(x)
I ('h
2
dx)
'R
ch
with (2.30)
" s m a l l " and
0.
2
W e use i n f a c t approximations o f
o b t a i n e d v i a a n u m e r i c a l i n t e g r a t i o n p r o c e.d u r e .
W e have t o m e n t i o n t h a t t h e o p t i m a l c h o i c e f o r
is still
Fh
a n open q u e s t i o n , n e v e r t h e l e s s t h e r e e x i s t s o m e e m p i r i c a l methods t h a t w e d o n ' t d i s c u s s h e r e .
Numerical r e s u l t s
ed w i t h t h e p i e c e w i s e l i n e a r a p p r o x i m a t i o n i n Sec. Eemark
(k=l)
obtairr
a r e given
2.6.
2.14:
T h e r e g u l a r i z a t i o n by
€
(or i t s
l ( h @ ) + 1 2 dx
v a r i a n t s ) m a y be s e e n a s a t e c h n i q u e f o r i n t r o d u c i n g s o m e pseudo-viscosity
B
-
i n t h e problem under c o n s i d e r a t i o n .
A method u s i n g ( 2 . 7 ) .
L e t u s d e s c r i b e t h i s method f o r t h e c o n t i n u o u s p r o b l e m , ~~
first:
We assume t h a t t h e e n t r o p y c o n d i t i o n i s f o r m u l a t e d by ( 2 . 7 ) . Let
M(x)
b e a s u f f i c i e n t l y s m o o t h u p p e r bound o f
e s t i m a t e d or g u e s s e d ) .
a weak f o r m u l a t i o n o f
R e p l a c i n g ( 2 . 7 ) by
(2.31)
is
hb,
(M i s
190
GLOWINSKI and 0. P I R O N N E A U
R.
-
(2.32)
1
I
I
[email protected] dx
n
M ( x ) c p dx
V
cp E 8 + ( n ) .
R
i t i s v e r y convenient t o i n t r o d u c e t h e
M
I n s t e a d of u s i n g
S
s o l u t i o n ( d e f i n e d modulo a c o n s t a n t ) o f
Ago = M
over
0,
= g
over
r,
(2.33)
then
i""'o.vv
(2.33) has t h e f ollow ing v a r i a t i o n a l formulation dx =
-
(
M(x)cp dx +
R
f
r
v cp t H 1 ( n ) ,
mdr
(2.34) I t f o l l o w s from ( 2 . 3 4 )
-i
(2.35)
W e observe t h a t
that
(2.32)
v ( O - @ O ) ' ~ c p dx
5
a (@-a0) = an
0.
c a n a l s o be w r i t t e n
v cp E
0
rs+(n).
F o r t h e d i s c r e t e problem t h e o b v i o u s s t r a t e g y seems t o be t h e following: F i r s t we a p p r o x i m a t e
vaoh-vvh
dx =
Q0
-
\
R
where
Mh
by
Go,,
solution of
Mh(x)vhdx +
f ghcphdr r
i s a c o n v e n i e n t a p p r o x i m a t i o n of
M.
+ e p , c vh,
191
LEAST SQUARE SOLUTION O F N O N L I N E A R PROBLEMS
( 2 . 3 5 ) by
F i n a l l y we a p p r o x i m a t e
1
-
(2.38)
n
'
v(@h-@o,)'vqh dx
I n f a c t t h e above " o b v i o u s ' l
qh
s t r a t e g y f o r t h e d i s c r e t e problem
h a s t o be m o d i f i e d f o r t h e f o l l o w i n g r e a s o n s : 1)
V:h
k=l,
If
can be g e n e r a t e d by t h e c a n o n i c a l b a s i s
2)
B u t i t i s n o t the c a s e f o r
Voh.
f u n c t i o n s of
k=l
C o m p u t a t i o n s u s i n g ( 2 . 3 8 ) and done w i t h t h a t t h e approximation o f
close t o
k=2. have shown
t h e s o l u t i o n i s n o t good
t h e p o i n t s a t which s o n i c l i n e s t o u c h
r.
To overcome t h e s e d i f f i c u l t i e s we may p r o c e e d a s f o l l o w s : If
k=l:
ed from 1
C,
Let
be t h e s e t o f
Nh,
t o
where
Vh,
c a n o n i c a l b a s i s of
Nh
t h e v e r t i c e s of
= dim(Vh).
Sh,number-
Wh
Let
be t h e
i.e.
(2.39) with
(2.40)
wi
E Vh,
W e observe t h a t
po - .s -i t i v e Z of
e
w
VL
i
2
of
Wi(Pj)
0
onver
Vh
v P j E c.,
= 6 1. J.
R,
+k
i,
i s g e n e r a t e d by
and t h a t
Bh.
the
Then i n s t e a d
t o formulate t h e d i s c r e t e entropy c o n d i t i o n
using (2.38)
one t a k e s
(2.41)
-
V(@h-@oh)*vp dh x
which i s e q u i v a l e n t t o t h e
s e t of
the
'
' Nh
ph
Vi
following l i n e a r
192
R.
GLOWINSKI a n d 0.
T-'IRO"EALl
i z c i u a 1i t y c o n s t r a i n t s
f
( 2 .I42 ) 5,
If
v(ql-cboh)
-1
-vwi
V
dx < 0
i = 1,
...,N h .
i s u s e d , i n s t e a d o€ ( 2 . 4 2 ) x e h a v e I
-
(2.'43)
v(~tl-@oll)'v\\idx
0
4
V
...,N h .
i = 1,
62 Computations done with
and u s i n g ( 2 . 4 2 )
k=l
have p r o d u c e d
good r e s u l t s . k=2: The t e c h n i q u e s
I.f ~
~~
t o be u s e d a r e m o r e c o m p l i c a t e d a n d
a r e d e t a i l e d i n G l o ~ i n s k i - P i r o n n c a u[ 1 1 ,
-Remark ..2.15:
(It holds f o r
If
k=1,2).
some D i r i c h l e t
r,
b o u n d a r y c o n d i t i o n s a r c ' p r e s c r i h e d somewhere o v e r
the
p o s i t i v e cone used t o de€irie t h e d i s c r e t e entropy c o n d i t i o n w i l l b e r e l a t e d t o t h e subspace of
Vtl
c o n s i s t i n g of
those
f u n c t i o n s v a n i s h i n g a t t h e boundary nodes corresponding t o the
(discrete) Dirichlet condition.
2.26:
Remark
(or
Mh)
may b e a n u n e a s y
computations. p a r t A,
Houever,
t~e c h n i q_u_e i s u s e d
e s p e c i a l l y for p r o f i l e
l i k e i n t h e r e g u l a r i z a t i o n method o f
t o handle
2.'+./+ A p p r o x.i m~-~ ation of ~
~~
(2.'+2) ( o r
if
Part A,
(2.17),
(2.18),
Xh
then
Xh
is
the d i s c r e t e entropy conditions.
2.4.3 e i t h e r
one u s e s t h e r e g u l a r i z a t i o n methods
or
~_
a penalty
(2.43)).
one u s e s t h e methods d e s c r i b e d i n S ec.
= Vh
if
X.
take i n t o account
e s s e n t i a l l y determined b y If
task,
some e m p i r i c a l r u l e s h a v e b e e n o b t a i n e d ,
If we d o n ' t
M
The o p t i m a l c h o i c e f o r the hound f u n c t i o n
of
Sec.
xh =
2.4.3,
i s d e f i n e d by t h e l i n e a r i n e q u a l i t y constraints
LEAST S Q U A R E SOLUTIOV O F NON L I N E A R PKOI)LE,PIS
(2.42)
2.5
( o r (2.43))
ICe-rative
~-
2.5.1
so
i f w e proceed
lution of
l i k e i n Sec.
193
2.4.3,
P a r t B.
t h e a p p~r o x i m a t e p r~. oblems.
-~
~~~
Preliminary statements.
G eneralitie - s.
We s h a l l r e s t r i c t o u r a t t e n t i o n t o t h e f o l l o w i n g s i t u a t i o n : *
R
i s hounded anti w e o n l y have Neumann boiindary
conditions.
~
We assume a l s o t h a t K u t t a - J o u k o w s k y required
i s not
since t h e i r treatment
conditions are not s p e c i f i c of
transonic
flows. '
We d o n ' t
take
i n t o account
*
We s u p p o s e t h a t
a
= 1
in
(2.17), ( 2 . 1 8 ) .
(2.14)
and t h a t an e n t r o p y
c o n d i t i o n i s f o r m u l a t e d u s i n g t h e method o f Part B,
with the formulation (2.43)
Sec.
2.4.3,
(control's constraints).
The c a s e where t h e r e g u l a r i z a t i o n method o f P a r t A i s u s e d i s s i m p l e r , and i n f a c t may b e s e e n a s a p a r t i c u l a r c a s e o f other ( a t least i f '
one u s e s p e n a l t y t o h a n d l e
the
(2.43)).
We f i n a l l y assume t h a t we u s e p i e c e w i s e l i n e a r f i n i t e elements
(k=l).
S i n c e we d o n ' t closed
take i n t o account
convex s u b s e t of
Vh
(2.17), ( 2 . 1 8 ) ,
d e f i n e d by
t h e a p p r o x i m a t e c o n t r o l problem
(2.43).
xh
i s the
Therefore
i s a non l i n e a r ( a n d non ~
c o n v e x ) programming p r o b l e m i n w h i c h t h e i n d e p e n d e n t v a r i a b l e
is
qh.
method.
~~
We s h a l l s o l v e t h i s p r o b l e m by a c o n j u g a t e g r a d i e n t ~~~
T o take i n t o account t h e l i n e a r c o n s t r a i n t s ( 2 . 4 3 )
we h a v e u s e d a p e n a l t y method
( f o r t h e s o l u t i o n by p e n a l t y -
d u a l i t y m e t h o d s v i a ____ augmented ._______ Lagrangians see GlowinskiP i r o n n e a u [ 11 )
.
194
R.
2.5.2
GLOWINSKI and 0 . PIRONNEAU
Penalty approximation of the discrete ____ problem. ~
~
We use the notation o f Sec. 2.4.3.
(Uh9Vh)h
-- _
vh
2
L (n)-scalar product: ~-
the following _ approximate (2.44)
Let us define over
Nh
c
3 i=l
miuh(Pi)vh(Pi),
pi E
ch,
tf
i,
where
i :f
m
-
6. =
measure union of
(nil, ni TE
=
a. 0
1'
6, such that Pi is a vertex of T;
the corresponding norm is denoted by
\*\,.
Then we define
It follows therefore from (2.45) that (2.43) can also be written
which is equivalent to
where in
cph
(2.48)'
(cp,)'
does not denote the positive part of
but the approximation of it defined by
LEAST SQUARE
In
(2.48)s
@Il
195
SOLUTION O F N O N LINEAR PROULEPIS
Eh
is a function of
t h-r~ ou _ g_h_ ( 2 . 2 3 )
and
.
r
i s a p o s i t i v e p a r a m e t er. Then w e c a n p r o v e t h e f o l l o w i n g :
- -~ P rop-osition
2 .1
-
If
t h e approximate c o n t r o l problem ~~
Min
(2.50) has a s o l u t i o n ,
Jh(Fh)
'h
I-1'
t h e n t h e s o l u t i o n o f - t h e p e n a l i z e d problem
(2.51) converges
Remark 2.17: ______
If
Ooh
= 0
notation) a variant of i n Sec.
2.4.3,
in
~~
the r e g u l a r i z e d f u n c t i o n a l i n t r o d u c e d
(2.51) i s a f i n i t e dimensional c o n t r o l
constraint.
We h a v e s o l v e d t h i s p r o b l e m
the n - on l i n e a r c o n j u g a t e
_g _r _ a d i e n t method
( s e e P o l a k [ 1, C h . 2 ,
more e f f e c t i v e ,
f o r our problem,
pp.
53-551)
w h i c h seems
than the F l e_ tche r.- R e e v e s _
The s c a l a r p r o d u c t u s e d i n t h i s a l g o r i t h m i s t h e
s c a l a r product important
+m.
(2.48) w e recover ( w i t h d i f f e r e~nt
u s i n g t h e P o l a k - R i b i & ~ ev e r s i o n o f
version.
+
P a r t A.
The p e n a l i z e d problem problem without
r
(2.50)
toasolutioLo-f
i n d u c e d by
H1(n)
over
s t e p i n t h e s o l u t i o n of
algorithm i s t h e computation of
(2.51)
Vh.
Therefore a very
by the ab o v e
the gradient
r _J h_
d5 h
~ (5,)
; this
196
R.
point i s discussed
GLOWINSKI and 0.
PIRONNEAU
i n t h e f o l l o w i n g s e c t i - o n 2.5.4.
A d e s c r i p t i o n o f t h e m o s t i m p o r t a n t non l i n e a r c o n j u g a t e g r a d i e n t a l g o r i t h m s i s g i v e n i n Sec.
2.5.4
Computation of
-- J h r
4.
( F ~ ) .
d5 h Owing t o t h e p r a c t i c a l i m p o r t a n c e of
.
i t s c o m p u t a t i o n w i t h some d e t a i l s :
Jhr
"h
we s h a l l d i s c u s s
We have
I
We c a n e a s i l y compute t h e l a s t t e r m of
(2.53). of
(2.23)
About
t h e right-hand
s i d e of
t h e s e c o n d t e r m , w e o b t a i n by d i f f e r e n t i a t i o n
197
LEAST SQUARE S O L U T I O N O F N O N L I N E A R PROBLEMS
Since
(2.57) !
t h e second term o f with t h e l a s t
Jhr ( 5 h )
' 6!3
t e r m of
~
= 0
t h e r i g h t hand s i d e o f
( 2 . 5 3 ) we o b t a i n
h'
Remark 2.18: I f -~ ..-.
a
( 2 . 5 6 ) i s e a s i l y computed and by a d d i t i o n
i n s t e a d of t a k i n g
U.=l
i n (2.24),
one t a k e s
w i l l r e q u i r e the use o f
t h e n t h e c o m p u t a t i o n of
an adjoint s t a t e equ ation. -~
2.5.5
Computational c o n s i d e r a t i o n s . __
When s o l v i n g ( 2 . 5 1 )
by t h e non l i n e a r c o n j u g a t e g r a d i e n t
method d i s c u s s e d a b o v e , we h a v e t o s o l v e a t e a c h i t e r a t i o n the state equation
(2.23).
S i n c e t h e b i l i n e a r form o c c u r i n g
i n (2.23)
is positive de€inite
p o i n t of
6
( o n c e t h e v a l u e of'
h a s b e e n p r e s c r i b e d ) we c a n u s e t o s o l v e
e i t h e r i-t e r a t i v e m e t h o d s l i k e c o n j u g a t e g r a d_--.I ient, I_
S.S.O.R.
etc...
A x e l s s o n [ 13
qjh
,
(cf.
e.g.,
Young [ 11 )
P o l a k [ 1 1 , Concus-Golub
i n one
(2.23)
S.O.R., __
[l]
,
o r d i r e c t m e-. thods l i k e G a u s s ' s o r
Choleskyls.
2.6
Numerical e x p e r i me n t s .
The f o l l o w i n g n u m e r i c a l t e s t s by Mrs. B r i s t e a u a n d MM.
( i n two d i m e n s i o n s ) were done
P e r i a u x , P o i r i e r on I B M 3 7 0 - 1 6 8 ,
198
R.
G L O W I N S K I and 0.
PIRONNEAU
o t h e r t e s t s and r e s u l t s a r e r e p o r t e d
C 11 ,
21
,
i n Glowinski-Pironneau
~ l o w i n s ~ < i - ~ e r i a u x - ~ i r o n n[c11 a u,[: 21
Concerning 3-dimensional
problems,
.
some c a l c u l a t i o n s h a v e been
d o n e on t h e e f f e c t s o f o b l i q u e wind on t h e f l o w a t t h e i n t a k e (see Fig.
2.4)
of
-
a j e t e n g i n e a t l o w mach number ; t h e s e
c a l c u l a t i o n s w i l l be r e p o r t e d
elsewhere.
Aircraft
Velocity
\
J e t engine
F i g u r e 2.4 The f o l l o w i n g r e s u l t s c o n c e r n
-
A Naca 0012 p r o f i l e ,
A Kornls p r o f i l e , An O n e r a ' s p r o f i l e ,
A two-piece
2.6.1
airfoil.
Naca 0012 p r o f i l e .
The t e s t d e s c r i b e d b e l o w c o r r e s p o n d s i n f i n i t y of
.85
(M,
= .85)
t o a mach number a t
and i s w i t h o u t _ i n~c_i _ d -e n c e . ~
We
199
L E A S T SQUARE SOLUTION OF NON LINEAR PROCLEMS
have used __piecewise linear finite elements ing to a triangulation %h
correspond-
(k=l)
of 600 triangles and 341 nodes.
The entropy condition (2.12) was required via regularization. ~
~~
The convergence was obtained in 80 iterations for a CPU time of
4
mn.
O n Figure 2.5 we have shown the mach numbers at the midpoint
of the edges, on the skin of the airfoil, extrapolated from their values at the centroid o f the triangles. ed by Bauer-Garabedian-Korn's
Results obtain-
code are shown also.
We observe
the good agreement between the results, except for the amplitude of the shock jump,
This seems to be related to the
fact that in the B-G-K code we have used the discretization -. finite difference - scheme is achieved by a Eon conservative -
and that the small disturbance model is used.
(The schemes
we use are fully conservative). O n Figure 2 . 6 we show the results obtained for the same problem, using piecewise quadratic finite elements 2.6.2
(k=2).
Korn's profile.
The Korn's profile, which is not symmetric (see Fig. 2.7) has the property of giving a shock free flow if
for s o m e specific incidence angle).
Mm = 0.75
(and
On Figure 2.8 we have
shown the mach numbers on the skin of the airfoil;
we observe
that the computed flow is also shock free. 2.6.3
ONERA profile.
Wind tunnel measurements have been made by ONERA (Office Nationale d'Etudes et de Recherches AGrospatiales) for the
R.
200
GLOWINSKI and 0 . P I R O N N E A U
profile shown o n Fig. 2.9.
O n Fig.
2.9 is also indicated the
mach distribution on the skin o f the profile.
O n Fig. 2.10
the computed results are shown and the agreement with the measured results is very good, including the shock location and amplitude.
The agreement is less good after the shock
since the physical flow is n o longer potential.
We observe
that the computed shock is sharper than the measured shock. This fact is due, may b e , to the fact the actual fluid is slightly viscous, but very likely the main reason of this difference is the fact that the velocity fluctuations at the intrance o f the wind-tunnel induce fluctuations of the shock wave location; therefore the measurement taking same average value of the velocity has a smoothing effect,
2.6.4
Two-piece airfoil. ____~
The geometry shown on Fig. 2.11 corresponds to a &)ractical situation.
Kutta-Joukowsky conditions have to be imposed o n
both profiles.
T h e position o f the front piece makes the
funnel slightly convergent. incidence and (k=l)
Mm = . 5 5 .
The main airfoil is at 5 O of
Lagrangian elements of order 1
were used with 2936 triangles and 1 5 5 5 nodes.
Regularization i s used and convergence is obtained in iterations for a C P U time o f 1 6 mn. guessed by measuring
3.5
X
from
lo3
and at
IlV(@-C)Il
n=50
2;
L it is 2.
T h e precision can be at
n=O
its value is
O n each triangle it varies
in the subsonic region to
sonic region.
50
lom2
i n the super-
li
LEAST SQUARE S O L U T I O N OF N O N L I N E A R PROBLEMS
i
0
20 1
202
0
*a
R. G L O W I N S K I and 0 . P I R O N N E A U
W
m
N
F
LEAST SQUARE SOLUTION OF NON L I N E A R I’R0DLJ)MS
203
204
a U I
R.
h
x
G L O W I N S K I and 0 . P I R O N N E A U
h
h
4
LEAST SQUARE SOLUTION O F N O N L I N E A R PROBLEMS
'I;
205
R.
206
G L O W I N S K I and 0 .
PIRONNEAU
M
X 0.5
Figure 2.10
1.
L E A S T SQUARF: S O L U T I O N O F NON L I N E A R PROBLEMS
R. GLOWINSKI and 0. PIRONNEAU
208
3.
Numerical solutions of Navier-Stokes equations. ___ __
3.1
-
Generalities.
Or$entati&n.
In this Section 3 we should like to briefly introduce the reader to some aspects o f the numerical analysis o f the NavierStokes equations.
This is a very large subject in constant
evolution; for information prior to 1 9 7 6 see Temarn [l] and the bibliography therein.
In the following subsections we s h a l l see how the non linear least square methodology can be applied for solving the continuous Navier-Stokes equations in the -~ pressure-velocity ~~
formulation.
Similar methods can also be used for the stream
function-vorticity __ formulation; these methods will bedescribed in a forthcoming report. In practice the above problems have to-be approximated;
we
have used some mixed finite element methods which are described in the report mentioned above.
However we shall only
consider the continuous case whose formalism is simpler.
3.2
The - _N _a v i e r - S t o k e s - e - q u a t i o n s .
Let
0
be a domain o f
iRN
(N = 2 , 3
in practice) and
we consider the
stationary Navier-Stokes equations
(3.2)
uIT. = g
where
(
g'n dr = 0 ,
Jl-
r
= 30;
2 01'
LLAST S Q U A R E S U 1 , I J T I O N O F N O N L I K E A R I'ROB1,CMS
F o r e x i s t e n c e anti u r i i q u c i i e s s see, e.g.,
7.3
Lions
I31 ,
(3.I)-( 3 . 3 )
tlieorems concerning
Tctnam [ 11
A l e a s t s q u.a.r e f o r m u l a t i o n
,
L a d y z e n s k a y a [ 11
.
o f~-t h e Y a v i e r - S t o k e s
equatio _ n_s .
It can b c
(3.4)
Miri VEU
where i n
v 5 :6
i
-
I v(y-v) I
dx,
'R
(3.4)
and where
i s a function of
v
through
(3.5)
W e h a v e o b t a i n e d a c-ontrol
problcrn i n w h i c h t h e s t a t e e q u a t i o n ~
~
~~
i s a Stokcs problem.
If t h e o r i g i n a l p r o b l e m
(3.1)-(3.3) has a s o l u t i o n then the
l e a s t s q u a r e p r o b l e m has a s o l u t i o n f o r w h i c h t h e c o s t function takes
the value
0.
Remayk 3 .-1-.: O t h e r f o r m u l a t i o n s c a n b e u s e d ; f o r i n s t a n c e we can r e q u i r e f o r
etc..
.
v
the condition
v-v
= 0,
o r t a k e as c o s t
R. GLOWINSKI and 0. PIRONNEAU
2 10
3.4
Calculation o f t h e gradient of the cost function. ___ - ._ ~. ~
~
~~
In practical applications, an essential step, in view of using descent or conjugate gradient methods, is the calculation o f the d e r i v a t i v e , o f the cost function. Let
J
J(v) = ; 2
where
r,
the functional defined by
y
is a function of
v
lv(y-v)l
via ( 3 . 5 ) .
2
dx,
Then we f o r m a l l y
have =
V(y-v).V&(y-v)dx
(3.6) = -V
where
6y
[
v(y-v)*06v
dx +
V
i
~ ( y - v *V6ydx )
W~V€(H’,(~))~
R
is solution o f
I
+
-VA6y
v6p = - ( & v . v ) v
-
(v-V)bv,
At that stage it is convenient to introduce the following S t o k e s equation (which is the adjoint state equation i n o u r
problem)
-van + m (3.8)
y,v E
(H;(anN
on
a,
= 0,
nlr
van = If
= -v~(y-v)
0.
(H1(n))N,
then (3.8) has a unique solution in
x (L2(W/R).
Taking variational formulations of ( 3 . 7 ) , (3.8) we have ~
LEAST SQUARE
where
Wo
Taking
z
SOLUTTON
= { z
E
= 7
i n ( 3 . 9 ) and
I n t e g r a t i n g by- p a r t hand s i d e o f
7 - 2 = 03
(Hi(n))N,
(i.e.
z
= 6y
211
L I N E A R PROBLEMS
O F NON
*
i n (3.10)
we o b t a i n
using Green's formula) the right-
( 3 . 1 1 ) c a n a l s o be w r i t t e n
I t f o l l o w s t h e n from ( 3 . 6 ) , ( 3 . 1 1 ) , ( 3 . 1 2 ) t h a t w 6 v E (H:(R))~
which c a n also be w r i t t e i i
(3.14) (
where and
(J'(V),~V)=(-VA(Y-V)
- ,*)
(H:(n))N.
(3.15)
3.5
+
(V*V)n
+
(V.V)q
-
(V~)q,q)
d e n o t e s t h e d u a l i t y p a i r i n g between
(H-l(n))N
Therefore
J' (V) =
-vA(V-y)
+
(V*V)
n
-k
(v'V:Y
-
(vV)n
-
Numerical e x p e r i m e n t s .
I n p r a c t i c e we have a p p r o x i m a t e d t h e N a v i e r - S t o k e s by a mixed
f i n i t e e l e m e n t method.
equations
Then u s i n g a d i s c r e t e
212
R.
v a r i a n t of
GLOWINSKI and 0 . P I R O N N E A U
the l e a s t square formulation
( 3 . ) + ) , ( 3 . 5 ) we have
minimized t h e c o s t f u n c t i o n by a c o n j u g a t e g r a d i e n t a l g o r i t h m _ _ ~ . ~ . _ ( o f R i b i & r e - P o l a k t s t y p e ) . I n t h e conjugate g r a d i e n t algorithm 1 N
we used as a s c a l a r P r o d u c t a d i s c r e t e f o r m of t h e ( ~ ~ ( 0 ) ) s c a l a r product.
The n u m e r i c a l t e s t p r e s e n t e d h e r e i s t h e t w o - d i m e n s i o n a l , time d e p e n d e n t f l o w p a s t a d i s k , o f a v i s c o u s i n c o m p r e s s i b l e f l u i d a t a Reynolds number o f 2 0 0 . f o r t h e time d i s c r e t i z a t i o n ,
U s i n g a n i_ m - p l i c i t scheme .__ ~
t h e numerical i n t e g r a t i o n i s
r e d u c e d t o a s e q u e n c e of s t a t i o n a r y a p p l y t h e t e c h n i q u e s . d e s c r i b e d above.
problems
t o which we
T h e r e f o r e a t each time
s t e p we h a v e t o s o l v e a l e a s t s q u a r e p r o b l e m ; however s i n c e
t o compute t h e f l o w
at
(n+l)At
we s t a r t o u r i t e r a t i v e
procedure u s i n g t h e informations a t of
nAt,
t h e convergence
the conjugate gradient algorithm i s f a i r l y f a s t .
S i n c e t h e domain o f t h e f l o w i s unbounded we have imbedded t h e d i s k i n a much l a r g e r d i s k on which a u n i f o r m boundary velocity i s prescribed. l y perturbated Fig.
3.1
S t a r t i n g from a u n i f o r m f l o w s l i g h t -
t o 3.4
show t h e e v o l u t i o n o f t h e f l o w
d u r i n g some t i m e i n t e r v a l , and p a r t i c u l a r l y t h e b e h a v i o u r of t h e f l o w behind t h e o b s t a c l e .
LEAST SQUARE
SOLUTION O F NOK L I N E A R PROBLEMS
F i g . 3.1
CYCLE DE TEMPS REYNOLDS
10
200.
214
R.
G L O W I N S K I and 0. P I R O N N E A U
Fig. 3 . 2
CYCLE DE TEMPS REY NClLDS
20
200.
LEAST S Q U A R E S O L U T I O N O F NON L I N E A R I’ROULEIIS
Fig. 3 . 3
CYCLE DE TEMPS REYNOLDS
40 200.
216
R.
G L O W I N S K I and 0 .
PIRONhTAU
Fig. 3 . 4
CYCLE DE TEMPS REYNOLDS
GO 200.
LEAST SQUARE S O L U T I O N OF NON L I N E A R PRODLEMS
4.
On c o n j u g a t e g r a d i e n t m e t h o d s .
An e s s e n t i a l t o o l i n t h e s o l u t i o n o f non l i n e a r p r o b l e m s by n j-u.g_ a t. e_ l e a s t s q u a r e methods i s t h e c_o_ - g r a d i e n t method t o b e d e s c r i b e d below:
N
be a s u f f i c i e n t l y smooth f u n c t i o n ,
Let
f: R
R
vf
af be {a fK ... , ,axN]
4
let
i t s g r a d i e n t , and l e t S b e a N x N 1 s y m m e t r i c , p o s i t i v e- d e_ f-i n i t e m a t r i x . Then t o s o l v e t h e m i n i =
~
~
r n i z a t i o n problem
I
x
Find
E RN
such t h a t
(4.1)
we c a n u s e t h e t w o f o l l o w i n g c o n j u g a t e g r a d i e n t a l g o r i t h m s t h e c o n v e r g e n c e o f which is s t u d i e d i n P o l a k [I]
(forS=I).
(Fletcher-Reeves) :
F i r s t method
(4.2)
xo
(4.3)
go =
s-l
(4.4)
w
-
-
0
E I R ~ given, Vf("),
0
= g .
Then a s s u m i n g t h a t
-
xn
-
wn
and
a r e known we compute
-
x n+1
by
(4.5)
X n + l= I
where problem
p,
-n
-p,y9
n
i s t h e s o l u t i o n f o r t h e one-dimensional minimization
R.
218
GLOWINSKI a n d 0 . PTRONNEAU
Then w e compute
-gn+'
(4.7)
gn+l =
(4.8)
wn+1 = -gn+1 + x n w n
wn+ 1
and
s-1
by
Vf(?n+l),
where
(4.9) i n (4.9),
denotes t h e u s u a l inner-product
( a , . )
Second method
of
N
R
.
(Polak-Ribisre).
T h i s method i s l i k e t h e p r e v i o u s o n e e x c e p t t h a t
(4.9) i s
r e p l a c e d by
-
( sgn+l, _g"+l-gn)
(4.10)
Remark 4.1:
If
E RN
where
and where
A
NxN
is a
symmetric p o s i t i v e
d e f i n i t e m a t r i x , t h e n t h e two a b o v e a l g o r i t h m s a r e i n f a c t
similar,
s i n c e i t c a n be proved t h a t i n t h i s c a s e
-
( S g n + l , g") Remark 4 . 2 :
A t e a c h i t e r a t i o n of
= 0.
t h e two a b o v e a l o r i t h m s w e
have t o s o l v e a one d i m e n s i o n a l problem t o o b t a i n
pn.
For
t h e s o l u t i o n o f t h i s k i n d o f one d i m e n s i o n a l problem s e e , e.g.,
P o l a k [: 1 3
Remark 4.3: round-off
,
Brent [
13 .
The above a l g o r i t h m s a r e f a i r l y s e n s i t i v e t o
e r r o r s ; h e n c e d o u b l e p r e c i s io n may be r e q u i r e d o n
some c o m p u t e r s .
M o r e o v e r i t may b e c o n v e n i e n t t o r e s t a r t
LEAST SQUARE SOLUTION OF MON LINEAR PROBLEMS
periodically with
wn =
-g
n
.
219
Actually Powell [l] has recently
introduced a new type of restarting procedure which increases the efficiency of the algorithm, but which requires more storage.
The Powell's analysis shows also that the Polak-
RibiBre's algorithm may b e more convenient than the FletcherReeves' algorithm for some problems.
In fact f o r the type of
problems we are solving a similar observation has been done, from our computational experiments. Remark 4.4: F o r the problems described i n this paper a convenient choice f o r the matrix form of
S
is to take a discrete
-A.
5 . Conclusions, The methods presented in this paper may appear as very heuristical, however when convexity o r monotonicity are lacking they seem one o f the few possibilities for computing the solutions o f some difficult non linear problems, taking into account their special structure. Applications to non linear problems i n elasticity are to be done. Acknowledgments: Many thanks are due to Mrs. Bri,steau and MM. Perrier, Periaux, Poirier for their help at various stages of this report.
R. GLOWINSKI and 0. PIRONNEAU
220
Bibliography Axelsson, 0.
-
[l]
A class of iterative methods for finite
element equations. Ccmp.
Math. . Applied _ ~ _ .. Mech. Eng.,
Bauer, F., Garabedian, P., Korn, D.
~~
-
[l]
Supercritical _____ wing
sections, . -~ Lecture Notes in Economics and Math.
.
Systems, Vol. 66, Springer-Verlag, 1972. Bauer, F., Garabedian, P., Korn, D., Jameson, A. [l]
Supercritical wing sections (11), Lecture Notes
in Economics and Math. Systems, Vol. 108, SpringerVerlag, 1975. Begis, D., Glowinski, R.-[l] 616ments finis
Application de la m6thode des la r6solution d'un problbme de
domaine optimal. Vol. 2, NO 2, Brent, R.
-
Applied Math. and Optimization,
(1975), pp. 130-169.
[l] Algorithms for minimization without derivat~
~~~
~
ives, __ Prentice Hall, 1973. Cea, J., Geymonat, G.
-
[l] Une methode de lin6arisation via
l'optimisation.
Instituto d i Alta Mat. Symp. Mat.,
Vol. X, (1972), Bologna, pp. Ciarlet, P.G., Raviart, P.A.
-
431-451.
111 Interpolation theory over
curved element with applications to finite elements methods. Comp. Meth. Applied Mech. Eng., Vol. 1, (1972) , PP. 217-249.
[2] A mixed finite element method for the biharmonic equation, i n Mathematical aspects of the
LEAST SQUARE SOLUTION OF NON LINEAR PROBLEMS
finite elements in partial di-fferen&ial C. de B o o r Ed., Acad. Press, Concus, P., Golub, G.H.
-
221
equations,
1974, pp. 125-145.
[l] A generalized conjugate gradient
method for non symmetric systems of linear equations, in Computing Methods in Applied Sciences and Engineering, R. Glowinski and J.L. Lions Ed., -~ Lecturca Notes i n Ecsiiomics and Math. Systems,
VOI. 134, Springer-Verlag, 1976, pp. 56-65. Glowinski, R.
-
[ 13 Numerical methods f o r non linear
variatipnal problems, Lecture Notes at the Tata Institute, Bombay, to appear in 1977-1978. [2] Approximations externes par 616ments finis
d'ordre un et deux d u probl&me de Dirichlet pour
A 2 , in Topics i n Numerical Analysis , J. J.H. Miller ed., Acad. Press, (1973)) pp. 123-171. Glowinski, R., Marroco, A .
-
[l] Sur lfapproximation par
616ments finis dfordre un, et la r6solution par p6nalisation-dualit6 d'une classe de probl;?mes de Dirichlet non linGaires, Revue Franqaise d'dutoma___ ____ tique, d'Informatique et de Recherche Opgr-ationnelle, e
9
ann&e, aojt 1975, R-2, pp.
41-76.
Glowinski, R., Periaux, J., Pironneau, 0.
-
[l] Transonic
flow simulation by the finite element method via optimal control.
Proceedings of the Second Inter-
national Symposium on Finite Elements i n Flow Problems, Santa Margherita, Italy, June 1976,
PP. 249-259.
R. GLOWINSKI and 0 . PIRONNEAU
222
[2] Use o f optimal control theory for the numerical
simulation of transonic flow by the method of finite lements.
Proceedings o f the
5th Conference _____
on Numerical Methods i n Fluid Dynamics.
Enschede,
Netherlands, July 1976. Glowinski, R., Pironneau, 0.
-
[l] On the computation of
transonic flows, Proceed+-ngs o f the first FrancoJ_po&Ese
Co-lloquium on Functional and Numerical
__ Analysis, Tokyo, Kyoto, September 1976. [2] Calcul d'ecoulements transoniques par des
methodes dlel6ments finis et de controle optimal, in Computing Methods in Applied Sciences and E n- -g i n e e r i x ,
R. Glowinski and J.L.
Lions Ed.,
Lectures Notes i n Economics and Math. Systems, vol. 134, Springer-Verlag,
1976, pp. 276-296.
[ 3 ] Numerical methods f o r the first biharmonic equation and for the two-dimensional- Stokes problem, C.S. Report, University of Stanford, Computer Science Dpt., Hewitt, B.L.;
1977.
Illingworth, C.R.;
Lock, R.C.;
Mc Donnel, J.H.; Richards, C.; Walkden, F . ;
Mangler, K.W.; Ed.
[l] Computational method2 and problems i n aeronau-
tical fluid Acad. Press, 1976. ~ - dynamics, -___-Jameson, A.
- [l]
Transonic flow calculations, in Proceedings
of Conference o n Computational Fluid Dynamics, - Von .
Karman Institute, March 1976, Brussels (Belgium). [2] Three dimensional flows around airfoils with
LEAST SQUARE SOLUTION OF NON LINEAR PROBLEMS
223
shoks , in Computing Methods _in ApplA-e-d-Sci-e-Eces_and Engineering, Part 2 , R . -~ _ I
Glowinski and J.L. Lions
Ed., Lecture Notes in Comp. Sciences, V o l . Springer-Verlag,
131
11,
(1974), pp. 185-212.
Iterative solution of transonic flows over
airfoils and wings, including flows at Mach 1,
Comm. Pure Appl. Math., ~ _ _Vol. 27, (1974), pp.283-309.
[4] Numerical solution of non linear partial differential equations of mixed type, in Numerical Solution of Equations, 111, Synspade - Partial Differential - __. ~~~
~~
1975, B. ___
~
Hubbard Ed., Acad. Press, (1976),
PP. 275-320. Landau, L.; Lifchitz, F.
-
[l] M6canique des - Fluides, Editions ~
MIR, MOSCOU, 1953. Ladyzenskaka, O . A .
-
[l] The ___. mathematical theory of viscous
incompressibA-e fluids, Gordon Breach, New York, 1963. Lions, J.L.
-
[l] Gontrole optimal des _ systemes _ ~ gouvern6s par ~
des 6quations aux d6rivGes partielles.
Dunod, 1968.
[ 2 ] Some aspects of the optimal control of
distributed systems, Regional Conference Series in Applied Math., Siam Publication NO 6, 1972.
[ 3 ] Quelques m6thodes de rdsolution des problhmes aux limites non lineaires, Dunod, 1969. Morawetz, C.S.
-
[l] O n the non existence of continuous
transonic flows past profiles, C o m m . Pure Applied Math., V o l .
9, 1956, pp. 45-68.
224
R. GLOWINSKI and 0. PIRONNEAU
Murman, E.M.;
Cole, J . D .
-
[11 Calculation of plane steady
transonic flows, AIAA Journal, V o l .
9, (1971),
pp. 114-121. Polak, E.
-
[l]
Computational methods in optimization, Acad.
Press, 1971. Powell, M . J . D .
-
[l] Restart procedures f o r the conjugate
gradient method, Math. Programming 12,
(1977),
pp. 241-254. Strang, G.; Fix, G.
-
[l] &n--analysis of the finite element -
method, Prentice Hall, 1973. Temam, R.
-
[l] Navier-Stokes Equations, North-Holland,
1977 Young, D . M .
- I: 11
Iterative s-ol~~ix-of__large linear systems,
Acad. Press, 1971.
LEA.I;T SQUARE SOLUTION O F NON L I N E A R PROBLEMS
I N F L U I D DYNAMICS
R.
GLOWINSKI* and
Universitb de P a r i s Iria-Laboria,
0.
PIRONNEAU
6 and
Iria-Laboria, 78150 le Chesnay
France
France
1. I n t r o d u c t i o n . Many problems i n Mathematical P h y s i c s a r e of
the
following forms, e i t h e r A(u)u = f,
+ B(u)
Au
In
(l.l), A ( v )
i s in (1.2)
v
i s for
=
f.
g i v e n a l_i _ n e_a_ r . o p e r a t o r and
B
a non l i n___ e a r operator.
This suggests obviously the following algorithms: F o r s o l v i n g (1.1)
(1.3)
uo
then f o r
n
2
0
we o b t a i n
u n+l
given from
u
n
by
I
A(un)un+' = f ,
For s o l v i n g ( 1 . 2 )
we may use an a l g o r i t h m s i m i l a r t o
(1.5)
(1.4)1 i s r e p l a c e d by
*
except t h a t
(1.3)-
~.
V i s i t i n g P r o f e s s o r IJF'RJ,
R i o d e J a n e i r o , July-August
1977.
172
R. GLOWINSKI and 0. P I R O N N E A U
There are some applications in which (l.3), (resp. ( 1 . 3 ) , ( 1 . 4 ) 2 , (resp. (1.2)).
-
(1.5))
( 1 . 4 ) 1’ (1.5)
converges to a solution o f (1.1)
Unfortunately these are important applications
mathematically modelled by equations like (1.1) or (1.2)
-
f o r which the above algorithms can not be used, therefore one
has to use more sophisticated methods.
In some cases (cf., e.g., Glowinski-Marocco Ell) these problems may be reduced to a saddle-point problem by introduc-. ~
ing an artificia-1 constraint and the corresponding Lagran* multiplier. - __ -..~
In some other situation such an approach is not suitable (by lack o f monotonicity usually) and it may be convenient to use aleast square method.
A typical least square formulation o f problem (1.1) is Min J(v)
(1.6)
V
where
(1.7) In
(1.7),
]I*]]
denotes a suitable norm (for more details
about this last point, see Remark 1.2 below) and
J(v)
the
solution o f the linear problem
We obtain similarly a least square formulation of (1.2) by replacing (1.8)1 by
(1*8)2
Ay = f-B(v).
LEAST SQUARE S O L U T I O N O F NON L I N E A R PROBLEMS
Remark _ _ 1~ .1:
I t a p p e a r s c l e a r l y from t h e above r e l a t i o n s t h a t
(1.6), (l.7), s t r u c t u r e of
(l.7),
( 1 . 8 ) l and ( 1 . 6 ) ,
( l . 8 ) 2 have t h e
p t i m a l c o n t r o l problem i n which t h e c ontrol an o ._ ~~
variable
(i.e.
t h e independent v a r i a b l e i n t h e problem) i s
the s t a t e variable i s
y
and t h e c o-s t-f-u n-c.t_ i -o. n is -
v,
J.
T h i s p a r t i c u l a r s t r u c t u r e l e a d s v e r y n a t u r a l l y t o compute t h e derivative o f
111 , [ 2 ] )
J
v i a an a d j o i n t s t a t e e q u a t i o n ( s e e Lions
f o r t h e o p t i m a l c o n t r o l o f s y s t e m s m o d e l l e d by
p a r t i a l d i f f e r e n t i a l equations). T h i s k i n d o f methodology h a s b e e n a d v o c a t e d by Cea-Geymonat
[l] f o r s o l v i n g non monotone, non l i n e a r t w o - p o i n t b o u n d a r y v a l u e problems l i k e
1
(1.9)
+
-u”
cp(u) = f
over
10,1[
u ( 0 ) = u ( 1 ) = 0.
I n B e g i s - G l o w i n s k i [1] a s i m i l a r a p p r o a c h h a s b e e n u s e d f o r s o l v i n g f r e e boundary problems
l i k e those r e l o t e d t o flows o f
l i q u i d s i n porous media. I n a n u n p u b l i s h e d work Derviauc-Glowinski h a v e a l s o u s e d t h i s t y p e o f methods t o compute t h e maximal s o l u t i o n of
I
(1.10)
where
>
-u
= Xe
U
over
10,1[
”
u(0) = u ( 1 ) = 0
0.
In t h e above r e f e r e n c e s t h e o p t i m i z a t i o n i s a c h i e v e d u s i n g e s s e n t i a l l y -~ descent o r steepest descent algorithms. ~
~~
~~~
More
~~~~
r e c e n t l y G l o w i n s k i and P i r o n n e a u o b s e r v e d t h a t t h e u s e of a conjugate gradient algorithm,
t o achieve the minimization,
c a n d r a s t i c a l l y improve t h e s p e e d o f c o n v e r g e n c e ,
specially
17 4
R. GLOWINSKI and 0. PIRONNEAU
in the solution of some hard problems like those to be considered in the following sections.
J
Remark 1.2: The choice of a convenient cost function
is a
very important matter if one wishes to use the above rnethodology.
Let us mention some possibilities we have found very
attractive since they simplify the adjoint state equation and lead to good convergence properties.
In the case of the problem (1.1) assume that for A(v)
E e(V,V’)
dual space. if
A(v)
where
(
where
is an Hilbert space and
V
Assume also thak
A(v)
v
given, V’
is self-adjoint. .
is _ strongly - a natural choice for _ _ _ V-elliptic
,*)
denotes the cu-ality -~_ pairing ~ _ . _o _ f _V’
Similarly i f in (1.2),
A
and
its Then
J
is
V.
is strongly V-elliptic it is
convenient to use
In the following sections we shall see how the above principles can be applied to the numerical solution of two important problems i n Fluid Mechanics, namely: *
The transonic flow of an isentropic, compressible perfect ~
fluid. ’ Navier-Stokes equations for incompressible fluids.
Numerical results will be presented in order to illustrate the possibilities of the above methods.
To conclude this report some indication will be given on the minimization of non quadratic functionals by conjugate gradient algorithms.
LEAST S Q U A R E S O L U T I O N OF NON L I N E A R PRODLEMS
175
2. Numerical solution-of transonic flow problems. . __ ~
2.1
Generalities on transonic flows.
The theoretical and numerical studies of transonic flows for -_ perfect fluids have always been very important,
But these
problems have become even more important these last years in relation with the design and development of large subsonic aircrafts. These transonic problem are very difficult, theoretically and mumerically, for the following reasons: They are n o n linear. Shocks may exist in the flow. An Entropy Condition is needed, in a way o r another, to avoid non physical solutions. They are mixed problems i n the sense that they are el_l.ipt-icin the Eubsgsic_ part of the flow, LyEerbolic i n the supersonic part of the flow. From a theoretical point of view let us mention the work of Morawetz [l].
At the present moment the numerical methods
which are the more commonly used have originated from Murman-Cole [l]
and we shall mention B a u e r - G a r a b e d i a n - K o r n r l ] ,
Bauer-Garabedian-Korn,Jameson
11
,
Jameson [ 11 ,[ 21 ,[ 31 ,[ 41
and the bibliographies there i n (see also Hewitt-Illingworth and Co-editors [ 13 ) . These above numerical methods use the key idea of Murman and Cole which consists to use a finite difference scheme, centered i n the subsonic-part of the flow, kackward (in the direction of the flow) in the supersonic part.
The switching
1 76
R.
GLOWINSRI a n d 0 .
between t h e s e two-schemes
PIRONNEAU
i s a u t o m a t i c a l l y done v i a a
t r u n c a t i o n o p e r a t o r only a c t i v e i n the supersanic p a r t of flow ( s e e Jameson,
loc.
cit.,
f o r more d e t a i l s ) .
We s h a l l d e s c r i b e h e r e a d i f f e r e n t a p p r o a c h ,
which seems v e r y
c o n v e n i e n t f o r comput i n g f l o w s p a s t p r o f i l e s ,
or f l o w s i n n o z z l e s .
infinity,
the
subsonic a t
I n o u r method t h e t r a n s o n i c
f l o w problem i s f o r m u l a t e d a s a non l i n e a r l e a s t s q u a r e problem,
control w h i c h i s i n t u r n i n t e r p r e t e d a s a n o p t i m a l~-
p r o b l -e m .
Then f i n i t e e l e m e n t s a r e u s e d t o a p p r o x i m a t e t h e Since t h e e n t r o p y~-c o n d i t i o n i s f o r m u l a t e d by .____
above problems.
~
a l i_ n e_ ar _ ine q u a l i t y c~-o n s t . ra i n t , a c o n v e n i e n t method t o h a n d l e .
i t i s t o u s e p e n a l t y a n d / o r d u a l i t y methods ( u s i n g a n augmented l a g r a n g i a n i f p e n a l t y and d u a l i t y a r e co mb in ed ) . Duality techniques a r e not
c o n s i d e r e d i n t h i s r e p o r t and w e
r e f e r t o Glowinski-Pironneau [ d e s c r i p t i o n of
2.2 2.2.1
11,
Glowinski
[13 for a
them.
M a t h e m a t i c a l m o d e l s f o r t r a n s__ o n __ i c f.l_ o ws. . .. B..-a s i c a s s u m p t i o n s .
We a s s u m e t h a t t h e f l u i d u n d e r c o n s i d e r a t i o n i s p e r f e c t a n d compressible,
and t h a t t h e flow of
a n d i r r o t a t i o-n a l ( i . e . ,
potential).
-
such a f l u i d is isentropic These assumptions a r e n o t
t r u e i n g e n e r a l s i n c e t h r o u g h a s h o c k wave t h e r e i s a v a r i a t i o n of e n t r o p y , and a n i r r o t a t i o n a l f l o w becomes rotational
(i.e.
t h e v a l i d i t y of i n t h e c a s e of
t h e r e i s c r e a t i o n o f v o_ r t_ i c i_ ty).
Therefore
t h e model t o f o l l o w i s assumed t o b e c o r r e c t "weak s h o cks". -.
177
LEAST SQUARE SOLUTION O F N O N LINEAR PROBLEMS
I n t h e c a s e o f a f l o w s p a s t a s- h a_ r p~ p_r o_f ._ i l e w e s h a l l suppose t h a t t h e r e i s no wake b e h i n d
2.2.2
E--q u a -~ t i o n s of
the tr a- _ i l ~i n g edge.
t h e flow.
~
n
Let
be t h e domain o f
r
t h e f l o w and
i t s boundary;
then
we s h a l l assume t h a t t h e f l o w i s m o d e l l e d by
-
(2.1)
v-((1 -
,
1
, 2
y-1
where i n ( 2 . 1 ) .
- 0 -
C,
-
y
the flow v e l o c i t y ,
i s t h e c.r i t i c a l v e l o c i t y , ~~
i s t h e r. __ a t i o of
W e have *
09
i s the flow p otential, _____
(y = 1 - 4 f o r
specific heats _ _ I
4.
t o add t o ( 2 . 1 )
Boundary c o n d i t i o n s ( o f D i r i c h l e t a n d / o r
Neumann t y p e , f o r
instance), *
-
A K u t t a - J o u k o v-k y
condition i n the case of
a l i f t i n g body ( s e e L a n d a u - L i f c h i t z
the flow past
[ l , Sec.
46]),
An e n t r o p y c o n d i t i o_ n . i n o r d e r t o e l i m i n a t e t h e non p h y s i c a l .s o l u t i o n s of
(1.1); t h i s p o i n t
w i l l be d i s c u s s e d i n S e c .
2.2.3.
:
Remark 2 . 1 :
I t c a n h a p p e n t h a t on some p a r t of
0
have
and
ness; -
t o be g i v e n s _ i m_ u l t a n e o u s- __ ly t o i n s u r e unique-
i t is, f o r i n s t a n c e ,
of F i g .
2.1 i f
~~
t h e c a s e of
'4
g i v e n on
t h e divergent nozzle
the velocity a t the entrance i s supersonic.
T y p i c a l boundary c o n d i t i o n s a r e 5-n
t h e boundary,
r4, rl, r 2 .
If
@
g i v e n on
r l , Fg
and
the velocity a t the entrance
178 (i.e.
R.
I-,)
GLOWINSKI and 0. P I R O N N E A U
is subsonic we require ___ less boundary conditions.
Figure 2.1
Figure 2.2 Remark 2.2:
In the case o f t h e f l o w past a multibody profile
(like i n Fig. 2.2)
each
body requires a Kutta-Jourkowsky
condition.
2.2.3
Formulation of the entropy condition.
It follows f r o m Landau-Lifchitz [ l , ch.91 condition c a n be formulated as follows:
that the entropy
1 79
LEAST SQUARE SOLUTION OF NON LINEAR PROBLEMS
1 In the direction
I
-!
(2.2)
o f the flow, One- cannot have a - ___ ~~
subsonic-supersonic transition through
a_;t?ck.
F o r a one-dimensional . . flow, (2.2) implies ~
*<
2 d
(2.3)
dx
3-
2
+m,
is a measure bounded from_ above; _ ~ _ _weak (and more dx rigorous) formulation of (2.3) are: i.e.
-~
There exists a constant ~~
~
M,
such that-either -
where
(We recall that fJ
(R) =
{cp E Cm ( f i )
,
cp
has a compact support in
n] . )
In the case o f a t ~ oo r three dimensional flow we shall suppose that (2.2) can be formulated by
(2.9)
190
R.
GLOWINSKI a n d 0 . P I R O N N E A U
T h e n u m e r i c a l r e s u l t s o b t a i n e d f o r t__ wo and _ t r e_e -_ d i_ m e~n s_ i o_ n a_ l flows u s i n g d i s c r e t e analogous of above formulation Remark 2 2 :
(2.7) seanto justify the
of t h e entropy condition.
The c o n d i t i o n ( 2 . 7 )
i s a c t u a l l y a p a r t i c u l a r case
O f
(w)+ E LPP)
(2.10) where,
assuming t h a t
is a m a s_ u_ re, _e_
( ~ c p ) + = positive p a r t of
(2.11) Then
Acp
(2.7)
corresponds t o
f r o m a computational point
p =
+m;
~cp.
another useful choice i s ,
of view,
p=2
i.c.
(2.12)
To c o n c l u d e t h i s S e c .
2.2
a t t h e moment,
l e t u s mention t h a t ,
t h e t r a n s o n i c f l o w problem i s s t i l l open from a t h e o r e t i c a l ~
~~
of view, even i n t h e simpler c a s e of t h e so-called
point
s m a l l d i s t u r b a n c_e _ c c l u a t i o n ~~
details).
( s e e J a m e s o n [ 11 ,[
41
f o r more
T h i s l a c k o f e x i s t e n c e and u n i q u e n e s s r e s u l t s h a s
v e r y l i k e l y i t s c o u n t e r p a r t i n t h e n u m e r i c a l methodology.
2.3
L - east
s q u a r e f o r m u l a t i on. ~
~~
R edu5;ion _ _t o~ a n o p t i m a l -~ ~~
~
~
~~
c o n t r o l problem. If
w e s u p p o s e t h a t t h e d e n s i t y o f t h e f l u i d i s o n e when t h e
velocity is zero,
then the coefficient
a s the density of
the
(2.13)
fluid.
of
V@
i n (2.1)
We s h a l l u s e t h e n o t a t i o n 1
lvcpl
2 Y - l
)
.
appm
181
LEAST SQUARE SOLUTION OF NON L I N E A R PROBLEXS
The i d e a o f t h e method t o f o l l o h i s t o d ecoup_ le t h e _ density ~ and t h e p o t e n t i a l those described
4
potential
0
and
-
via a _ l e_a s -~ t s_ q u_ a r_ e~ f o_ rm u_ l a.t_____ i o_ n like _ _
@
To d o s o we i n t r o d u c e a new
1.
i n Sec.
-
the c -o .n t r o l potential ~
by m i n i m i z i n g a c o n v e n i e n t
least
square c o s t function.
W e may u s e f o r i n s t a n c e t h e f o l l o w i n g f o r m u l a t i o n formulations s e e Glowinski-Pironneau pironneau
5
and t r y t o r e c o u p l e
(for other
[ 2 1 , Gowinski-Periaux-
[ 11 )
(2.14) where
(2.15) i
In
(2.14
,
a
'' c o n v e n i e n t l y ''
p ep(x)) = 0
y =
that for a i r , appears t h a t
via the state
p l u s boundary c o n d i t i o n s f o r
t h e parameter
a convex s e t Since
5
i s a f u n c t i o n of
in
i n the
is either cho s e n
.
implies
1
y+1 2
(--)
Y-1
t r a n s o n i c range
=
1
= (Fl) --
,,/6 = 2 . 4 5 , =
(say
X
and
1 Y+1 2
I Vcp(x) 1
i f and o n l y i f
1.4
or
0
I-.
on
d,
is
and
C,
it
lVCpl
1.5
c*)
w e have
(2.16)
o
< 6 s p(cp(x)) s 1
I t f o l l o w s from ( 2 . 1 6 )
that
a.e.
on
0.
( 2 . 1 5 ) i s an e l l i p t i c p r o b l e m f o r
a p p r o p r i a t e boundary c o n d i t i o n s . l i f t i n g b o d i e s , Kutta-Joukowsky
I n t h e case of flows p a s t conditions are a l s o required
physical t o h a v e , w i t h t h e o t h e r b o u n d a r y c o n d i t i o n s , a -__ s o l u t i o n o f problem ( 2 . 1 5 )
(modulo a c o n s t a n t i f
Neumann c o n d i t i o n s on t h e b o u n d a r y ) .
one o n l y h a s
R.
182
2.4:
Remark
If
GLOWINSKI and 0 . P I R O N N E A U
i n t h e o r i g i n a l problem,
0
r,
be s i m u l t a n e o u s l y p r e s c r i b e d on some p a r t of approach,
5
w i t h t h e two p o t e n t i a l
and
'@ have t o an
and
a,
t h e above
i s very
convenient s i n c e t h e boundary c o n d i t i o n s can be s p l i t t e d b e t ween
@
4.
and
However i f
boundary c o n d i t i o n s f o r t a k e a c c o u n t of
0
one w i s h e s t o u s e t h e same and
4
i t i s always p o s s i b l e t o
t h e e x t r a boundary c o n d i t i o n s
(assumed t o be
(2.14) a
o f D i r i c h l e t t y p e ) by a d d i n g t o t h e c o s t f u n c t i o n I
quantity proportional t o e i t h e r
( I
' rd
15-@,1
(9-0,1
2
br,
or
'd (or t o a l i n e a r c o m b i n a t i o n o f b o t h ) , where
2
i s the part
of
on which one r e q u i r e s
A s i m i l a r i d e a i s used i n Begis-Glowinski
@Ird
=
rd
Od.
[l] t o s o l v e some
f r e e boundary p r o b l e m s . Remark 2.5: variant
To s t a t e t h e e n t r o p y c__ o n d-i t_i_ o n ( 2 . 7 )
4
( 2 . 1 2 ) ) we h a v e t h e c h o i c e between
(or its
and
F
(actud-
l y we can a l s o u s e t h e s e two p o t e n t i a l s ~ i m y l t a n e o u ~ l y ) I. f one u s e s
E
( r e s p . @ ) we have a c o n t r o l c o-n s t r a_i_n t
(resp.
a s t a t e c-onstr-int). Remark 2 . 6 :
W e observe t h a t for t h e c l a s s o f flows under
c o n s i d e r a t i o n we h a v e
T h i s o b s e r v a t i o n s u g g e s t s t h a t t h e convex
X
( 2 . 1 4 ) h a s t o be t a k e n a s a convex s u b s e t of
occuring i n
Wl'"(S2).
We
o b s e r v e a l s o t h a t t o s t a y i n t h e t r a n s o n i c r a n g e i t may b e convenient t o intr oduce t h e following c o n s t r a i n t s
or
183
LEAST SQUARE S O L U T I O N OF NON LINEAR PROBLEMS
(2.18) Actually the computations we have done prove that for a physically well-posed transonic problem it is not necessary ~~
to introduce
Remark
~
(2.17), (2.18).
2.7: If the transonic problem has a solution and if
X
is "large enough" then the least square problem has a solution such that the cost function is equal to zero.
This last
property gives us means to check the quality of the computed solution.
2.4
Finite element approximations. __
We assume that
0c
_. Three-dimensional calculations will
2
!R
be reparted i n a forthcoming paper.
2.4.1
Generalities.
The above control problem will be approximated by finite elements.
Compared to finite differences, finite elements
give us the possibility o f handling problems posed on rather complicated geometrical domains.
Moreover the variational . .~
framework of finite element formulations is very appropriate to the problem under consideration and to the methodology we use to solve it.
It will be in particular fairly easy to
approximate the weak formulations of the entropy condition (2.7) (and also the alternate entropy condition (2.12)).
If
R
is unbounded, it is replaced by a bounded domain
still denoted by
-
as large as possible.
-
To approximate
the above continuous problem, we introduce a standard tri-
184
R. GLOWINSKI and 0 . PIRONNEAU
Zh
angulat ion elements).
of
61
(we can also use quadrilateral finite -
7
Then we approximate the functions
and
$I
by
piecewise polynomial functions belonging to the following subspace of
with
H1(n)
(and
Wl’”(n))
Pk = the space o f polynomials of degree
Sk.
Remark 2.8: In the case of a lifting body, to take into account the Kutta-Joukowsky condition, one usually introduces -~ .. (see Fig. 2.3) an arc
y
between the trailing edge o f the
body and the external boundary.
Thus arc
constant jump, a priori unknown of
a@
requires the continuity of
Y.
Since
@
-
an
@
(and
is discontinuous along
used any more, but introducing
.
.
0
and
y
supports a
<,
ae z) at y,
but one the crossing o f
Hl(61)
c a n not be
\
n = n - y one can use
A.
HL(n)
and define
Vh
over a triangulation of
Results of computations for such profiles are given in
Sec. 2.6;
however the method f o r t r e a t i n g t h e Kutta-Joukowsky
condition will not be treated in this report since they are not specific of transonic flows. In order to simplify the presentation o f our method we shall assume i n the sequel that one works directly over
61.
LEAST SQUARF S O L U T I O \
@rU O U
Figure
Remark 2 . 9 :
If
@
and
?
2.3
have t o s a t i s f y o n l y N e u m-a r i r i
b-o u n d a r y_ _c o n d i t -~ i o n s , t h e s u b s p a c e of ~
Vh
itself.
If
on t h e c o n t r a r )
i$
Vh
t o h e u s e d w i l l be
and/or
haxe to satisfi
F
D i r i c h l e t boundary c o n d i t i o n s over s o m e p a r t of s h a l l have t o use subspaccs o f Remark 2 . 1 0 : domain o f
R
Vh,
,
",
t h e n we
s t r i c t l y included _-___
We have t a c i t l y assumed t h a t 2
185
LIhFAR l'ROI3LEVS
R
in V
h'
i s a polygonal
or has b e e n a p p r o x i m a t e d b y s u c h a d o m a i n .
Hoizever i n t h e c a s e o f a curx-ed b o u n d a r y i t i s a l w a y s r > o s s i b l e t o use
( a t s o m e e x t r a computational c o s t ) curved f i n i t e
elements
t-11).
(see,
for i n s t a n c e , S t r a r i g - F i x [I], C i a r l e t - R a v i a r t
R. GLOWINSKI and 0. PIRONNEAU
185
2.4.2
o f the state equation . and . . o f_ _the cost Approximation __. ~~
function. T o simplify the presention, we shall assume that we only have
Neumann .boundary cond>t-ions, i.e.
aa@ n
(2.20)
= g
on
7
(it is the case for the important application of flows subsonic at infinity, past profiles). F
sections that
and
0
It follows from the above
will satisfy the same boundary ~
conditions.
We shall also assume that if
then the corresponding value of (2.21)
p g
p
dr =
g(x)
f 0, x E
r,
is known, and that 0.
'iThen the state equation (2.15) h a s the following variational
(2.22)
which is approximated by
r
i
4
(2.23) where @
and
P(ch)v@h'v"h
f 'I-
(pg)h q h
[
Oh vh' ( ~ g )is~ a Convenient approximation of @h
dr
pg.
qhevh,
Since
are only defined modulo a constant we shall
prescribe the values of some point of ed by
dx =
r.
@
and
cph
(and
5
and
5,)
at
The cost function in (2.14) is approximat-
187
LEAST SQUARE S Q L U T I O N O F PJON L I N E A R PROULEFIS
If one
R e m a r k 2.11:
(k=l),
uses
the piecewise linear approximation
then the integrals occuring i n ( 2 . 2 3 ) ,
to compute since property for
05,
p(ch)
(2.24)
piecewise constant implies the same (from (2.1'3)).
If
k=2
a ~ numerical _
integration procedure has to be used in ( 2 . 2 3 ) ,
a
are easy
_
and also in
= 1.
(2.24)
if
2.4.3
Approximation of .the entropy condition. -~ _. ~~
~
~~~
~
T o avoid non physical shocks (i.e. shocks for which the
thermodynamical entropy condition is n o t satisfied) we use the formulations o f Sec. 2 . 2 . 3 .
We still assume that w e o n l y
have Neumanii b o u n d a r y conditions.
We use the notation of the continuous oroblem. the condition ( 2 . 1 2 ) (2.14)
If
one uses
we may "regularize" the cost function
by adding to it, with
C
> 0, either
(2.25)
or (2.26)
E
[ '*
l ( h E ; ) + 1 2 dx
( o r a linear combination of both).
In
(2.25),
(2.26),
is
a "small" parameter.
We can make this approach more sophisticated by using instead of ( 2 . 2 5 ) ,
(2.26),
regularization functionals like
188
R.
GLOWINSKI and 0 . PIRONNEAU
0.
p o s s i b l y e q u a l t o z e r o o v e r some p a r t of
To u s e t h e above methodology f o r t h e a p p r o x i m a t e problem i t i s n e c e s s a r y t o have a n a p p r o x i m a t i o n o f W e may u s e a n a p p r o x i m a t i o n a p p r o x i m a t i o n of
C i a r l e t - R a v i a r t [ 21
,
[
(2.28)
73 ) :
i s s u f f i c i e n t l y smooth, t h e n from
we h a v e
Ai/lcpdx =
'n
(2
i
epdP
-
(2.29)
I
{R
an a p p r o x i m a t i o n of
E H1(n).
-
i
'R
Ah
$I,
0Bh'VCPh dx
A@
of
cph E Vh,
bL
I n (2.29),
t o define
the function
of
r,
g
of
gh
Remark 2.13: To o b t a i n
(2.20).
t h i s can be d o n e . from ( 2 . 2 9 ) we have t o s o l v e a
Ah@,
l-~ i n e a r s ys t e m t h e m a t r i x o f which i s _s_yrnmetric, p o s i t i v e ~~
sparse, but p o t d i a g o n a l
approximation of t h e o p e r a t o r approximate
is
t h e above i d e a i s s l i g h t l y more
complicated t o a p p l y , b u t
definite,
as
we a l s o h a v e D i r i c h l e t boundary c o n d i t i o n s
If
o v e r some p a r t
cp
E Vh*
We u s e t h e same method
Remark 2 . 1 2 :
ghCPhc
=
'h'hcphdx
AhOh
+
V w - v c p dx
F
U s i n g t h i s i d e a we d e f i n e a n a p p r o x i m a t i o n €allows
[21,
e q u a t i o n ( s e e Glowinski
-
Glowinski-Pironneau [
$
L e t u s assume t h a t Grgeq-for5qJ-a-
s u g g e s t e d by F i x e d -f i n i t e e_l e_m_e n t
t h e b i h a r m o n~ic
~
AS).
(resp.
I
n
4 0 cp
h h h
dx
I).
( t h i s m a t r i x i s an If
k=l
one c a n
u s i n g t h e two-dimensional,
189
LEAST SQUARE SOLUTION O F N O N L I N E A R I'RODLEMS
t____ r a p e-z o i d a l n u m e r i c a l i n t e g r a t i o - n m e t h o d .
D o i n g s o we o b t a i n
~
s o l v i n g a l i n e a r s y s t e m w i t h a _d_i a g o n a l m a t r i x .
by
ahQh
If
t h e a b o v e r e g u l a r i z a t i o n method i s t e c h n i c a l l y more
k=2
complicated t o use.
A
Once
h a s b e e n a p p r o x i m a t e d we add t o t h e c o s t f u n c t i o n
(see ( 2 . 2 4 ) ) the functional
2
I (Ahah)+]
('. 3')
dx
(rcsp'
I
Eh(x)
I ('h
2
dx)
'R
ch
with (2.30)
" s m a l l " and
0.
2
W e use i n f a c t approximations o f
o b t a i n e d v i a a n u m e r i c a l i n t e g r a t i o n p r o c e.d u r e .
W e have t o m e n t i o n t h a t t h e o p t i m a l c h o i c e f o r
is still
Fh
a n open q u e s t i o n , n e v e r t h e l e s s t h e r e e x i s t s o m e e m p i r i c a l methods t h a t w e d o n ' t d i s c u s s h e r e .
Numerical r e s u l t s
ed w i t h t h e p i e c e w i s e l i n e a r a p p r o x i m a t i o n i n Sec. Eemark
(k=l)
obtairr
a r e given
2.6.
2.14:
T h e r e g u l a r i z a t i o n by
€
(or i t s
l ( h @ ) + 1 2 dx
v a r i a n t s ) m a y be s e e n a s a t e c h n i q u e f o r i n t r o d u c i n g s o m e pseudo-viscosity
B
-
i n t h e problem under c o n s i d e r a t i o n .
A method u s i n g ( 2 . 7 ) .
L e t u s d e s c r i b e t h i s method f o r t h e c o n t i n u o u s p r o b l e m , ~~
first:
We assume t h a t t h e e n t r o p y c o n d i t i o n i s f o r m u l a t e d by ( 2 . 7 ) . Let
M(x)
b e a s u f f i c i e n t l y s m o o t h u p p e r bound o f
e s t i m a t e d or g u e s s e d ) .
a weak f o r m u l a t i o n o f
R e p l a c i n g ( 2 . 7 ) by
(2.31)
is
hb,
(M i s
190
GLOWINSKI and 0. P I R O N N E A U
R.
-
(2.32)
1
I
I
[email protected] dx
n
M ( x ) c p dx
V
cp E 8 + ( n ) .
R
i t i s v e r y convenient t o i n t r o d u c e t h e
M
I n s t e a d of u s i n g
S
s o l u t i o n ( d e f i n e d modulo a c o n s t a n t ) o f
Ago = M
over
0,
= g
over
r,
(2.33)
then
i""'o.vv
(2.33) has t h e f ollow ing v a r i a t i o n a l formulation dx =
-
(
M(x)cp dx +
R
f
r
v cp t H 1 ( n ) ,
mdr
(2.34) I t f o l l o w s from ( 2 . 3 4 )
-i
(2.35)
W e observe t h a t
that
(2.32)
v ( O - @ O ) ' ~ c p dx
5
a (@-a0) = an
0.
c a n a l s o be w r i t t e n
v cp E
0
rs+(n).
F o r t h e d i s c r e t e problem t h e o b v i o u s s t r a t e g y seems t o be t h e following: F i r s t we a p p r o x i m a t e
vaoh-vvh
dx =
Q0
-
\
R
where
Mh
by
Go,,
solution of
Mh(x)vhdx +
f ghcphdr r
i s a c o n v e n i e n t a p p r o x i m a t i o n of
M.
+ e p , c vh,
191
LEAST SQUARE SOLUTION O F N O N L I N E A R PROBLEMS
( 2 . 3 5 ) by
F i n a l l y we a p p r o x i m a t e
1
-
(2.38)
n
'
v(@h-@o,)'vqh dx
I n f a c t t h e above " o b v i o u s ' l
qh
s t r a t e g y f o r t h e d i s c r e t e problem
h a s t o be m o d i f i e d f o r t h e f o l l o w i n g r e a s o n s : 1)
V:h
k=l,
If
can be g e n e r a t e d by t h e c a n o n i c a l b a s i s
2)
B u t i t i s n o t the c a s e f o r
Voh.
f u n c t i o n s of
k=l
C o m p u t a t i o n s u s i n g ( 2 . 3 8 ) and done w i t h t h a t t h e approximation o f
close t o
k=2. have shown
t h e s o l u t i o n i s n o t good
t h e p o i n t s a t which s o n i c l i n e s t o u c h
r.
To overcome t h e s e d i f f i c u l t i e s we may p r o c e e d a s f o l l o w s : If
k=l:
ed from 1
C,
Let
be t h e s e t o f
Nh,
t o
where
Vh,
c a n o n i c a l b a s i s of
Nh
t h e v e r t i c e s of
= dim(Vh).
Sh,number-
Wh
Let
be t h e
i.e.
(2.39) with
(2.40)
wi
E Vh,
W e observe t h a t
po - .s -i t i v e Z of
e
w
VL
i
2
of
Wi(Pj)
0
onver
Vh
v P j E c.,
= 6 1. J.
R,
+k
i,
i s g e n e r a t e d by
and t h a t
Bh.
the
Then i n s t e a d
t o formulate t h e d i s c r e t e entropy c o n d i t i o n
using (2.38)
one t a k e s
(2.41)
-
V(@h-@oh)*vp dh x
which i s e q u i v a l e n t t o t h e
s e t of
the
'
' Nh
ph
Vi
following l i n e a r
192
R.
GLOWINSKI a n d 0.
T-'IRO"EALl
i z c i u a 1i t y c o n s t r a i n t s
f
( 2 .I42 ) 5,
If
v(ql-cboh)
-1
-vwi
V
dx < 0
i = 1,
...,N h .
i s u s e d , i n s t e a d o€ ( 2 . 4 2 ) x e h a v e I
-
(2.'43)
v(~tl-@oll)'v\\idx
0
4
V
...,N h .
i = 1,
62 Computations done with
and u s i n g ( 2 . 4 2 )
k=l
have p r o d u c e d
good r e s u l t s . k=2: The t e c h n i q u e s
I.f ~
~~
t o be u s e d a r e m o r e c o m p l i c a t e d a n d
a r e d e t a i l e d i n G l o ~ i n s k i - P i r o n n c a u[ 1 1 ,
-Remark ..2.15:
(It holds f o r
If
k=1,2).
some D i r i c h l e t
r,
b o u n d a r y c o n d i t i o n s a r c ' p r e s c r i h e d somewhere o v e r
the
p o s i t i v e cone used t o de€irie t h e d i s c r e t e entropy c o n d i t i o n w i l l b e r e l a t e d t o t h e subspace of
Vtl
c o n s i s t i n g of
those
f u n c t i o n s v a n i s h i n g a t t h e boundary nodes corresponding t o the
(discrete) Dirichlet condition.
2.26:
Remark
(or
Mh)
may b e a n u n e a s y
computations. p a r t A,
Houever,
t~e c h n i q_u_e i s u s e d
e s p e c i a l l y for p r o f i l e
l i k e i n t h e r e g u l a r i z a t i o n method o f
t o handle
2.'+./+ A p p r o x.i m~-~ ation of ~
~~
(2.'+2) ( o r
if
Part A,
(2.17),
(2.18),
Xh
then
Xh
is
the d i s c r e t e entropy conditions.
2.4.3 e i t h e r
one u s e s t h e r e g u l a r i z a t i o n methods
or
~_
a penalty
(2.43)).
one u s e s t h e methods d e s c r i b e d i n S ec.
= Vh
if
X.
take i n t o account
e s s e n t i a l l y determined b y If
task,
some e m p i r i c a l r u l e s h a v e b e e n o b t a i n e d ,
If we d o n ' t
M
The o p t i m a l c h o i c e f o r the hound f u n c t i o n
of
Sec.
xh =
2.4.3,
i s d e f i n e d by t h e l i n e a r i n e q u a l i t y constraints
LEAST S Q U A R E SOLUTIOV O F NON L I N E A R PKOI)LE,PIS
(2.42)
2.5
( o r (2.43))
ICe-rative
~-
2.5.1
so
i f w e proceed
lution of
l i k e i n Sec.
193
2.4.3,
P a r t B.
t h e a p p~r o x i m a t e p r~. oblems.
-~
~~~
Preliminary statements.
G eneralitie - s.
We s h a l l r e s t r i c t o u r a t t e n t i o n t o t h e f o l l o w i n g s i t u a t i o n : *
R
i s hounded anti w e o n l y have Neumann boiindary
conditions.
~
We assume a l s o t h a t K u t t a - J o u k o w s k y required
i s not
since t h e i r treatment
conditions are not s p e c i f i c of
transonic
flows. '
We d o n ' t
take
i n t o account
*
We s u p p o s e t h a t
a
= 1
in
(2.17), ( 2 . 1 8 ) .
(2.14)
and t h a t an e n t r o p y
c o n d i t i o n i s f o r m u l a t e d u s i n g t h e method o f Part B,
with the formulation (2.43)
Sec.
2.4.3,
(control's constraints).
The c a s e where t h e r e g u l a r i z a t i o n method o f P a r t A i s u s e d i s s i m p l e r , and i n f a c t may b e s e e n a s a p a r t i c u l a r c a s e o f other ( a t least i f '
one u s e s p e n a l t y t o h a n d l e
the
(2.43)).
We f i n a l l y assume t h a t we u s e p i e c e w i s e l i n e a r f i n i t e elements
(k=l).
S i n c e we d o n ' t closed
take i n t o account
convex s u b s e t of
Vh
(2.17), ( 2 . 1 8 ) ,
d e f i n e d by
t h e a p p r o x i m a t e c o n t r o l problem
(2.43).
xh
i s the
Therefore
i s a non l i n e a r ( a n d non ~
c o n v e x ) programming p r o b l e m i n w h i c h t h e i n d e p e n d e n t v a r i a b l e
is
qh.
method.
~~
We s h a l l s o l v e t h i s p r o b l e m by a c o n j u g a t e g r a d i e n t ~~~
T o take i n t o account t h e l i n e a r c o n s t r a i n t s ( 2 . 4 3 )
we h a v e u s e d a p e n a l t y method
( f o r t h e s o l u t i o n by p e n a l t y -
d u a l i t y m e t h o d s v i a ____ augmented ._______ Lagrangians see GlowinskiP i r o n n e a u [ 11 )
.
194
R.
2.5.2
GLOWINSKI and 0 . PIRONNEAU
Penalty approximation of the discrete ____ problem. ~
~
We use the notation o f Sec. 2.4.3.
(Uh9Vh)h
-- _
vh
2
L (n)-scalar product: ~-
the following _ approximate (2.44)
Let us define over
Nh
c
3 i=l
miuh(Pi)vh(Pi),
pi E
ch,
tf
i,
where
i :f
m
-
6. =
measure union of
(nil, ni TE
=
a. 0
1'
6, such that Pi is a vertex of T;
the corresponding norm is denoted by
\*\,.
Then we define
It follows therefore from (2.45) that (2.43) can also be written
which is equivalent to
where in
cph
(2.48)'
(cp,)'
does not denote the positive part of
but the approximation of it defined by
LEAST SQUARE
In
(2.48)s
@Il
195
SOLUTION O F N O N LINEAR PROULEPIS
Eh
is a function of
t h-r~ ou _ g_h_ ( 2 . 2 3 )
and
.
r
i s a p o s i t i v e p a r a m e t er. Then w e c a n p r o v e t h e f o l l o w i n g :
- -~ P rop-osition
2 .1
-
If
t h e approximate c o n t r o l problem ~~
Min
(2.50) has a s o l u t i o n ,
Jh(Fh)
'h
I-1'
t h e n t h e s o l u t i o n o f - t h e p e n a l i z e d problem
(2.51) converges
Remark 2.17: ______
If
Ooh
= 0
notation) a variant of i n Sec.
2.4.3,
in
~~
the r e g u l a r i z e d f u n c t i o n a l i n t r o d u c e d
(2.51) i s a f i n i t e dimensional c o n t r o l
constraint.
We h a v e s o l v e d t h i s p r o b l e m
the n - on l i n e a r c o n j u g a t e
_g _r _ a d i e n t method
( s e e P o l a k [ 1, C h . 2 ,
more e f f e c t i v e ,
f o r our problem,
pp.
53-551)
w h i c h seems
than the F l e_ tche r.- R e e v e s _
The s c a l a r p r o d u c t u s e d i n t h i s a l g o r i t h m i s t h e
s c a l a r product important
+m.
(2.48) w e recover ( w i t h d i f f e r e~nt
u s i n g t h e P o l a k - R i b i & ~ ev e r s i o n o f
version.
+
P a r t A.
The p e n a l i z e d problem problem without
r
(2.50)
toasolutioLo-f
i n d u c e d by
H1(n)
over
s t e p i n t h e s o l u t i o n of
algorithm i s t h e computation of
(2.51)
Vh.
Therefore a very
by the ab o v e
the gradient
r _J h_
d5 h
~ (5,)
; this
196
R.
point i s discussed
GLOWINSKI and 0.
PIRONNEAU
i n t h e f o l l o w i n g s e c t i - o n 2.5.4.
A d e s c r i p t i o n o f t h e m o s t i m p o r t a n t non l i n e a r c o n j u g a t e g r a d i e n t a l g o r i t h m s i s g i v e n i n Sec.
2.5.4
Computation of
-- J h r
4.
( F ~ ) .
d5 h Owing t o t h e p r a c t i c a l i m p o r t a n c e of
.
i t s c o m p u t a t i o n w i t h some d e t a i l s :
Jhr
"h
we s h a l l d i s c u s s
We have
I
We c a n e a s i l y compute t h e l a s t t e r m of
(2.53). of
(2.23)
About
t h e right-hand
s i d e of
t h e s e c o n d t e r m , w e o b t a i n by d i f f e r e n t i a t i o n
197
LEAST SQUARE S O L U T I O N O F N O N L I N E A R PROBLEMS
Since
(2.57) !
t h e second term o f with t h e l a s t
Jhr ( 5 h )
' 6!3
t e r m of
~
= 0
t h e r i g h t hand s i d e o f
( 2 . 5 3 ) we o b t a i n
h'
Remark 2.18: I f -~ ..-.
a
( 2 . 5 6 ) i s e a s i l y computed and by a d d i t i o n
i n s t e a d of t a k i n g
U.=l
i n (2.24),
one t a k e s
w i l l r e q u i r e the use o f
t h e n t h e c o m p u t a t i o n of
an adjoint s t a t e equ ation. -~
2.5.5
Computational c o n s i d e r a t i o n s . __
When s o l v i n g ( 2 . 5 1 )
by t h e non l i n e a r c o n j u g a t e g r a d i e n t
method d i s c u s s e d a b o v e , we h a v e t o s o l v e a t e a c h i t e r a t i o n the state equation
(2.23).
S i n c e t h e b i l i n e a r form o c c u r i n g
i n (2.23)
is positive de€inite
p o i n t of
6
( o n c e t h e v a l u e of'
h a s b e e n p r e s c r i b e d ) we c a n u s e t o s o l v e
e i t h e r i-t e r a t i v e m e t h o d s l i k e c o n j u g a t e g r a d_--.I ient, I_
S.S.O.R.
etc...
A x e l s s o n [ 13
qjh
,
(cf.
e.g.,
Young [ 11 )
P o l a k [ 1 1 , Concus-Golub
i n one
(2.23)
S.O.R., __
[l]
,
o r d i r e c t m e-. thods l i k e G a u s s ' s o r
Choleskyls.
2.6
Numerical e x p e r i me n t s .
The f o l l o w i n g n u m e r i c a l t e s t s by Mrs. B r i s t e a u a n d MM.
( i n two d i m e n s i o n s ) were done
P e r i a u x , P o i r i e r on I B M 3 7 0 - 1 6 8 ,
198
R.
G L O W I N S K I and 0.
PIRONNEAU
o t h e r t e s t s and r e s u l t s a r e r e p o r t e d
C 11 ,
21
,
i n Glowinski-Pironneau
~ l o w i n s ~ < i - ~ e r i a u x - ~ i r o n n[c11 a u,[: 21
Concerning 3-dimensional
problems,
.
some c a l c u l a t i o n s h a v e been
d o n e on t h e e f f e c t s o f o b l i q u e wind on t h e f l o w a t t h e i n t a k e (see Fig.
2.4)
of
-
a j e t e n g i n e a t l o w mach number ; t h e s e
c a l c u l a t i o n s w i l l be r e p o r t e d
elsewhere.
Aircraft
Velocity
\
J e t engine
F i g u r e 2.4 The f o l l o w i n g r e s u l t s c o n c e r n
-
A Naca 0012 p r o f i l e ,
A Kornls p r o f i l e , An O n e r a ' s p r o f i l e ,
A two-piece
2.6.1
airfoil.
Naca 0012 p r o f i l e .
The t e s t d e s c r i b e d b e l o w c o r r e s p o n d s i n f i n i t y of
.85
(M,
= .85)
t o a mach number a t
and i s w i t h o u t _ i n~c_i _ d -e n c e . ~
We
199
L E A S T SQUARE SOLUTION OF NON LINEAR PROCLEMS
have used __piecewise linear finite elements ing to a triangulation %h
correspond-
(k=l)
of 600 triangles and 341 nodes.
The entropy condition (2.12) was required via regularization. ~
~~
The convergence was obtained in 80 iterations for a CPU time of
4
mn.
O n Figure 2.5 we have shown the mach numbers at the midpoint
of the edges, on the skin of the airfoil, extrapolated from their values at the centroid o f the triangles. ed by Bauer-Garabedian-Korn's
Results obtain-
code are shown also.
We observe
the good agreement between the results, except for the amplitude of the shock jump,
This seems to be related to the
fact that in the B-G-K code we have used the discretization -. finite difference - scheme is achieved by a Eon conservative -
and that the small disturbance model is used.
(The schemes
we use are fully conservative). O n Figure 2 . 6 we show the results obtained for the same problem, using piecewise quadratic finite elements 2.6.2
(k=2).
Korn's profile.
The Korn's profile, which is not symmetric (see Fig. 2.7) has the property of giving a shock free flow if
for s o m e specific incidence angle).
Mm = 0.75
(and
On Figure 2.8 we have
shown the mach numbers on the skin of the airfoil;
we observe
that the computed flow is also shock free. 2.6.3
ONERA profile.
Wind tunnel measurements have been made by ONERA (Office Nationale d'Etudes et de Recherches AGrospatiales) for the
R.
200
GLOWINSKI and 0 . P I R O N N E A U
profile shown o n Fig. 2.9.
O n Fig.
2.9 is also indicated the
mach distribution on the skin o f the profile.
O n Fig. 2.10
the computed results are shown and the agreement with the measured results is very good, including the shock location and amplitude.
The agreement is less good after the shock
since the physical flow is n o longer potential.
We observe
that the computed shock is sharper than the measured shock. This fact is due, may b e , to the fact the actual fluid is slightly viscous, but very likely the main reason of this difference is the fact that the velocity fluctuations at the intrance o f the wind-tunnel induce fluctuations of the shock wave location; therefore the measurement taking same average value of the velocity has a smoothing effect,
2.6.4
Two-piece airfoil. ____~
The geometry shown on Fig. 2.11 corresponds to a &)ractical situation.
Kutta-Joukowsky conditions have to be imposed o n
both profiles.
T h e position o f the front piece makes the
funnel slightly convergent. incidence and (k=l)
Mm = . 5 5 .
The main airfoil is at 5 O of
Lagrangian elements of order 1
were used with 2936 triangles and 1 5 5 5 nodes.
Regularization i s used and convergence is obtained in iterations for a C P U time o f 1 6 mn. guessed by measuring
3.5
X
from
lo3
and at
IlV(@-C)Il
n=50
2;
L it is 2.
T h e precision can be at
n=O
its value is
O n each triangle it varies
in the subsonic region to
sonic region.
50
lom2
i n the super-
li
LEAST SQUARE S O L U T I O N OF N O N L I N E A R PROBLEMS
i
0
20 1
202
0
*a
R. G L O W I N S K I and 0 . P I R O N N E A U
W
m
N
F
LEAST SQUARE SOLUTION OF NON L I N E A R I’R0DLJ)MS
203
204
a U I
R.
h
x
G L O W I N S K I and 0 . P I R O N N E A U
h
h
4
LEAST SQUARE SOLUTION O F N O N L I N E A R PROBLEMS
'I;
205
R.
206
G L O W I N S K I and 0 .
PIRONNEAU
M
X 0.5
Figure 2.10
1.
L E A S T SQUARF: S O L U T I O N O F NON L I N E A R PROBLEMS
R. GLOWINSKI and 0. PIRONNEAU
208
3.
Numerical solutions of Navier-Stokes equations. ___ __
3.1
-
Generalities.
Or$entati&n.
In this Section 3 we should like to briefly introduce the reader to some aspects o f the numerical analysis o f the NavierStokes equations.
This is a very large subject in constant
evolution; for information prior to 1 9 7 6 see Temarn [l] and the bibliography therein.
In the following subsections we s h a l l see how the non linear least square methodology can be applied for solving the continuous Navier-Stokes equations in the -~ pressure-velocity ~~
formulation.
Similar methods can also be used for the stream
function-vorticity __ formulation; these methods will bedescribed in a forthcoming report. In practice the above problems have to-be approximated;
we
have used some mixed finite element methods which are described in the report mentioned above.
However we shall only
consider the continuous case whose formalism is simpler.
3.2
The - _N _a v i e r - S t o k e s - e - q u a t i o n s .
Let
0
be a domain o f
iRN
(N = 2 , 3
in practice) and
we consider the
stationary Navier-Stokes equations
(3.2)
uIT. = g
where
(
g'n dr = 0 ,
Jl-
r
= 30;
2 01'
LLAST S Q U A R E S U 1 , I J T I O N O F N O N L I K E A R I'ROB1,CMS
F o r e x i s t e n c e anti u r i i q u c i i e s s see, e.g.,
7.3
Lions
I31 ,
(3.I)-( 3 . 3 )
tlieorems concerning
Tctnam [ 11
A l e a s t s q u.a.r e f o r m u l a t i o n
,
L a d y z e n s k a y a [ 11
.
o f~-t h e Y a v i e r - S t o k e s
equatio _ n_s .
It can b c
(3.4)
Miri VEU
where i n
v 5 :6
i
-
I v(y-v) I
dx,
'R
(3.4)
and where
i s a function of
v
through
(3.5)
W e h a v e o b t a i n e d a c-ontrol
problcrn i n w h i c h t h e s t a t e e q u a t i o n ~
~
~~
i s a Stokcs problem.
If t h e o r i g i n a l p r o b l e m
(3.1)-(3.3) has a s o l u t i o n then the
l e a s t s q u a r e p r o b l e m has a s o l u t i o n f o r w h i c h t h e c o s t function takes
the value
0.
Remayk 3 .-1-.: O t h e r f o r m u l a t i o n s c a n b e u s e d ; f o r i n s t a n c e we can r e q u i r e f o r
etc..
.
v
the condition
v-v
= 0,
o r t a k e as c o s t
R. GLOWINSKI and 0. PIRONNEAU
2 10
3.4
Calculation o f t h e gradient of the cost function. ___ - ._ ~. ~
~
~~
In practical applications, an essential step, in view of using descent or conjugate gradient methods, is the calculation o f the d e r i v a t i v e , o f the cost function. Let
J
J(v) = ; 2
where
r,
the functional defined by
y
is a function of
v
lv(y-v)l
via ( 3 . 5 ) .
2
dx,
Then we f o r m a l l y
have =
V(y-v).V&(y-v)dx
(3.6) = -V
where
6y
[
v(y-v)*06v
dx +
V
i
~ ( y - v *V6ydx )
W~V€(H’,(~))~
R
is solution o f
I
+
-VA6y
v6p = - ( & v . v ) v
-
(v-V)bv,
At that stage it is convenient to introduce the following S t o k e s equation (which is the adjoint state equation i n o u r
problem)
-van + m (3.8)
y,v E
(H;(anN
on
a,
= 0,
nlr
van = If
= -v~(y-v)
0.
(H1(n))N,
then (3.8) has a unique solution in
x (L2(W/R).
Taking variational formulations of ( 3 . 7 ) , (3.8) we have ~
LEAST SQUARE
where
Wo
Taking
z
SOLUTTON
= { z
E
= 7
i n ( 3 . 9 ) and
I n t e g r a t i n g by- p a r t hand s i d e o f
7 - 2 = 03
(Hi(n))N,
(i.e.
z
= 6y
211
L I N E A R PROBLEMS
O F NON
*
i n (3.10)
we o b t a i n
using Green's formula) the right-
( 3 . 1 1 ) c a n a l s o be w r i t t e n
I t f o l l o w s t h e n from ( 3 . 6 ) , ( 3 . 1 1 ) , ( 3 . 1 2 ) t h a t w 6 v E (H:(R))~
which c a n also be w r i t t e i i
(3.14) (
where and
(J'(V),~V)=(-VA(Y-V)
- ,*)
(H:(n))N.
(3.15)
3.5
+
(V*V)n
+
(V.V)q
-
(V~)q,q)
d e n o t e s t h e d u a l i t y p a i r i n g between
(H-l(n))N
Therefore
J' (V) =
-vA(V-y)
+
(V*V)
n
-k
(v'V:Y
-
(vV)n
-
Numerical e x p e r i m e n t s .
I n p r a c t i c e we have a p p r o x i m a t e d t h e N a v i e r - S t o k e s by a mixed
f i n i t e e l e m e n t method.
equations
Then u s i n g a d i s c r e t e
212
R.
v a r i a n t of
GLOWINSKI and 0 . P I R O N N E A U
the l e a s t square formulation
( 3 . ) + ) , ( 3 . 5 ) we have
minimized t h e c o s t f u n c t i o n by a c o n j u g a t e g r a d i e n t a l g o r i t h m _ _ ~ . ~ . _ ( o f R i b i & r e - P o l a k t s t y p e ) . I n t h e conjugate g r a d i e n t algorithm 1 N
we used as a s c a l a r P r o d u c t a d i s c r e t e f o r m of t h e ( ~ ~ ( 0 ) ) s c a l a r product.
The n u m e r i c a l t e s t p r e s e n t e d h e r e i s t h e t w o - d i m e n s i o n a l , time d e p e n d e n t f l o w p a s t a d i s k , o f a v i s c o u s i n c o m p r e s s i b l e f l u i d a t a Reynolds number o f 2 0 0 . f o r t h e time d i s c r e t i z a t i o n ,
U s i n g a n i_ m - p l i c i t scheme .__ ~
t h e numerical i n t e g r a t i o n i s
r e d u c e d t o a s e q u e n c e of s t a t i o n a r y a p p l y t h e t e c h n i q u e s . d e s c r i b e d above.
problems
t o which we
T h e r e f o r e a t each time
s t e p we h a v e t o s o l v e a l e a s t s q u a r e p r o b l e m ; however s i n c e
t o compute t h e f l o w
at
(n+l)At
we s t a r t o u r i t e r a t i v e
procedure u s i n g t h e informations a t of
nAt,
t h e convergence
the conjugate gradient algorithm i s f a i r l y f a s t .
S i n c e t h e domain o f t h e f l o w i s unbounded we have imbedded t h e d i s k i n a much l a r g e r d i s k on which a u n i f o r m boundary velocity i s prescribed. l y perturbated Fig.
3.1
S t a r t i n g from a u n i f o r m f l o w s l i g h t -
t o 3.4
show t h e e v o l u t i o n o f t h e f l o w
d u r i n g some t i m e i n t e r v a l , and p a r t i c u l a r l y t h e b e h a v i o u r of t h e f l o w behind t h e o b s t a c l e .
LEAST SQUARE
SOLUTION O F NOK L I N E A R PROBLEMS
F i g . 3.1
CYCLE DE TEMPS REYNOLDS
10
200.
214
R.
G L O W I N S K I and 0. P I R O N N E A U
Fig. 3 . 2
CYCLE DE TEMPS REY NClLDS
20
200.
LEAST S Q U A R E S O L U T I O N O F NON L I N E A R I’ROULEIIS
Fig. 3 . 3
CYCLE DE TEMPS REYNOLDS
40 200.
216
R.
G L O W I N S K I and 0 .
PIRONhTAU
Fig. 3 . 4
CYCLE DE TEMPS REYNOLDS
GO 200.
LEAST SQUARE S O L U T I O N OF NON L I N E A R PRODLEMS
4.
On c o n j u g a t e g r a d i e n t m e t h o d s .
An e s s e n t i a l t o o l i n t h e s o l u t i o n o f non l i n e a r p r o b l e m s by n j-u.g_ a t. e_ l e a s t s q u a r e methods i s t h e c_o_ - g r a d i e n t method t o b e d e s c r i b e d below:
N
be a s u f f i c i e n t l y smooth f u n c t i o n ,
Let
f: R
R
vf
af be {a fK ... , ,axN]
4
let
i t s g r a d i e n t , and l e t S b e a N x N 1 s y m m e t r i c , p o s i t i v e- d e_ f-i n i t e m a t r i x . Then t o s o l v e t h e m i n i =
~
~
r n i z a t i o n problem
I
x
Find
E RN
such t h a t
(4.1)
we c a n u s e t h e t w o f o l l o w i n g c o n j u g a t e g r a d i e n t a l g o r i t h m s t h e c o n v e r g e n c e o f which is s t u d i e d i n P o l a k [I]
(forS=I).
(Fletcher-Reeves) :
F i r s t method
(4.2)
xo
(4.3)
go =
s-l
(4.4)
w
-
-
0
E I R ~ given, Vf("),
0
= g .
Then a s s u m i n g t h a t
-
xn
-
wn
and
a r e known we compute
-
x n+1
by
(4.5)
X n + l= I
where problem
p,
-n
-p,y9
n
i s t h e s o l u t i o n f o r t h e one-dimensional minimization
R.
218
GLOWINSKI a n d 0 . PTRONNEAU
Then w e compute
-gn+'
(4.7)
gn+l =
(4.8)
wn+1 = -gn+1 + x n w n
wn+ 1
and
s-1
by
Vf(?n+l),
where
(4.9) i n (4.9),
denotes t h e u s u a l inner-product
( a , . )
Second method
of
N
R
.
(Polak-Ribisre).
T h i s method i s l i k e t h e p r e v i o u s o n e e x c e p t t h a t
(4.9) i s
r e p l a c e d by
-
( sgn+l, _g"+l-gn)
(4.10)
Remark 4.1:
If
E RN
where
and where
A
NxN
is a
symmetric p o s i t i v e
d e f i n i t e m a t r i x , t h e n t h e two a b o v e a l g o r i t h m s a r e i n f a c t
similar,
s i n c e i t c a n be proved t h a t i n t h i s c a s e
-
( S g n + l , g") Remark 4 . 2 :
A t e a c h i t e r a t i o n of
= 0.
t h e two a b o v e a l o r i t h m s w e
have t o s o l v e a one d i m e n s i o n a l problem t o o b t a i n
pn.
For
t h e s o l u t i o n o f t h i s k i n d o f one d i m e n s i o n a l problem s e e , e.g.,
P o l a k [: 1 3
Remark 4.3: round-off
,
Brent [
13 .
The above a l g o r i t h m s a r e f a i r l y s e n s i t i v e t o
e r r o r s ; h e n c e d o u b l e p r e c i s io n may be r e q u i r e d o n
some c o m p u t e r s .
M o r e o v e r i t may b e c o n v e n i e n t t o r e s t a r t
LEAST SQUARE SOLUTION OF MON LINEAR PROBLEMS
periodically with
wn =
-g
n
.
219
Actually Powell [l] has recently
introduced a new type of restarting procedure which increases the efficiency of the algorithm, but which requires more storage.
The Powell's analysis shows also that the Polak-
RibiBre's algorithm may b e more convenient than the FletcherReeves' algorithm for some problems.
In fact f o r the type of
problems we are solving a similar observation has been done, from our computational experiments. Remark 4.4: F o r the problems described i n this paper a convenient choice f o r the matrix form of
S
is to take a discrete
-A.
5 . Conclusions, The methods presented in this paper may appear as very heuristical, however when convexity o r monotonicity are lacking they seem one o f the few possibilities for computing the solutions o f some difficult non linear problems, taking into account their special structure. Applications to non linear problems i n elasticity are to be done. Acknowledgments: Many thanks are due to Mrs. Bri,steau and MM. Perrier, Periaux, Poirier for their help at various stages of this report.
R. GLOWINSKI and 0. PIRONNEAU
220
Bibliography Axelsson, 0.
-
[l]
A class of iterative methods for finite
element equations. Ccmp.
Math. . Applied _ ~ _ .. Mech. Eng.,
Bauer, F., Garabedian, P., Korn, D.
~~
-
[l]
Supercritical _____ wing
sections, . -~ Lecture Notes in Economics and Math.
.
Systems, Vol. 66, Springer-Verlag, 1972. Bauer, F., Garabedian, P., Korn, D., Jameson, A. [l]
Supercritical wing sections (11), Lecture Notes
in Economics and Math. Systems, Vol. 108, SpringerVerlag, 1975. Begis, D., Glowinski, R.-[l] 616ments finis
Application de la m6thode des la r6solution d'un problbme de
domaine optimal. Vol. 2, NO 2, Brent, R.
-
Applied Math. and Optimization,
(1975), pp. 130-169.
[l] Algorithms for minimization without derivat~
~~~
~
ives, __ Prentice Hall, 1973. Cea, J., Geymonat, G.
-
[l] Une methode de lin6arisation via
l'optimisation.
Instituto d i Alta Mat. Symp. Mat.,
Vol. X, (1972), Bologna, pp. Ciarlet, P.G., Raviart, P.A.
-
431-451.
111 Interpolation theory over
curved element with applications to finite elements methods. Comp. Meth. Applied Mech. Eng., Vol. 1, (1972) , PP. 217-249.
[2] A mixed finite element method for the biharmonic equation, i n Mathematical aspects of the
LEAST SQUARE SOLUTION OF NON LINEAR PROBLEMS
finite elements in partial di-fferen&ial C. de B o o r Ed., Acad. Press, Concus, P., Golub, G.H.
-
221
equations,
1974, pp. 125-145.
[l] A generalized conjugate gradient
method for non symmetric systems of linear equations, in Computing Methods in Applied Sciences and Engineering, R. Glowinski and J.L. Lions Ed., -~ Lecturca Notes i n Ecsiiomics and Math. Systems,
VOI. 134, Springer-Verlag, 1976, pp. 56-65. Glowinski, R.
-
[ 13 Numerical methods f o r non linear
variatipnal problems, Lecture Notes at the Tata Institute, Bombay, to appear in 1977-1978. [2] Approximations externes par 616ments finis
d'ordre un et deux d u probl&me de Dirichlet pour
A 2 , in Topics i n Numerical Analysis , J. J.H. Miller ed., Acad. Press, (1973)) pp. 123-171. Glowinski, R., Marroco, A .
-
[l] Sur lfapproximation par
616ments finis dfordre un, et la r6solution par p6nalisation-dualit6 d'une classe de probl;?mes de Dirichlet non linGaires, Revue Franqaise d'dutoma___ ____ tique, d'Informatique et de Recherche Opgr-ationnelle, e
9
ann&e, aojt 1975, R-2, pp.
41-76.
Glowinski, R., Periaux, J., Pironneau, 0.
-
[l] Transonic
flow simulation by the finite element method via optimal control.
Proceedings of the Second Inter-
national Symposium on Finite Elements i n Flow Problems, Santa Margherita, Italy, June 1976,
PP. 249-259.
R. GLOWINSKI and 0 . PIRONNEAU
222
[2] Use o f optimal control theory for the numerical
simulation of transonic flow by the method of finite lements.
Proceedings o f the
5th Conference _____
on Numerical Methods i n Fluid Dynamics.
Enschede,
Netherlands, July 1976. Glowinski, R., Pironneau, 0.
-
[l] On the computation of
transonic flows, Proceed+-ngs o f the first FrancoJ_po&Ese
Co-lloquium on Functional and Numerical
__ Analysis, Tokyo, Kyoto, September 1976. [2] Calcul d'ecoulements transoniques par des
methodes dlel6ments finis et de controle optimal, in Computing Methods in Applied Sciences and E n- -g i n e e r i x ,
R. Glowinski and J.L.
Lions Ed.,
Lectures Notes i n Economics and Math. Systems, vol. 134, Springer-Verlag,
1976, pp. 276-296.
[ 3 ] Numerical methods f o r the first biharmonic equation and for the two-dimensional- Stokes problem, C.S. Report, University of Stanford, Computer Science Dpt., Hewitt, B.L.;
1977.
Illingworth, C.R.;
Lock, R.C.;
Mc Donnel, J.H.; Richards, C.; Walkden, F . ;
Mangler, K.W.; Ed.
[l] Computational method2 and problems i n aeronau-
tical fluid Acad. Press, 1976. ~ - dynamics, -___-Jameson, A.
- [l]
Transonic flow calculations, in Proceedings
of Conference o n Computational Fluid Dynamics, - Von .
Karman Institute, March 1976, Brussels (Belgium). [2] Three dimensional flows around airfoils with
LEAST SQUARE SOLUTION OF NON LINEAR PROBLEMS
223
shoks , in Computing Methods _in ApplA-e-d-Sci-e-Eces_and Engineering, Part 2 , R . -~ _ I
Glowinski and J.L. Lions
Ed., Lecture Notes in Comp. Sciences, V o l . Springer-Verlag,
131
11,
(1974), pp. 185-212.
Iterative solution of transonic flows over
airfoils and wings, including flows at Mach 1,
Comm. Pure Appl. Math., ~ _ _Vol. 27, (1974), pp.283-309.
[4] Numerical solution of non linear partial differential equations of mixed type, in Numerical Solution of Equations, 111, Synspade - Partial Differential - __. ~~~
~~
1975, B. ___
~
Hubbard Ed., Acad. Press, (1976),
PP. 275-320. Landau, L.; Lifchitz, F.
-
[l] M6canique des - Fluides, Editions ~
MIR, MOSCOU, 1953. Ladyzenskaka, O . A .
-
[l] The ___. mathematical theory of viscous
incompressibA-e fluids, Gordon Breach, New York, 1963. Lions, J.L.
-
[l] Gontrole optimal des _ systemes _ ~ gouvern6s par ~
des 6quations aux d6rivGes partielles.
Dunod, 1968.
[ 2 ] Some aspects of the optimal control of
distributed systems, Regional Conference Series in Applied Math., Siam Publication NO 6, 1972.
[ 3 ] Quelques m6thodes de rdsolution des problhmes aux limites non lineaires, Dunod, 1969. Morawetz, C.S.
-
[l] O n the non existence of continuous
transonic flows past profiles, C o m m . Pure Applied Math., V o l .
9, 1956, pp. 45-68.
224
R. GLOWINSKI and 0. PIRONNEAU
Murman, E.M.;
Cole, J . D .
-
[11 Calculation of plane steady
transonic flows, AIAA Journal, V o l .
9, (1971),
pp. 114-121. Polak, E.
-
[l]
Computational methods in optimization, Acad.
Press, 1971. Powell, M . J . D .
-
[l] Restart procedures f o r the conjugate
gradient method, Math. Programming 12,
(1977),
pp. 241-254. Strang, G.; Fix, G.
-
[l] &n--analysis of the finite element -
method, Prentice Hall, 1973. Temam, R.
-
[l] Navier-Stokes Equations, North-Holland,
1977 Young, D . M .
- I: 11
Iterative s-ol~~ix-of__large linear systems,
Acad. Press, 1971.