1. No category

Overview
Introduction
Theory
Implementation
Modifications
Results
Description of AdjointShapeOptimizationFoam
and how to implement new cost functions
Ulf Nilsson
A project work in the course
CFD with OpenSource software,
taught by
H˚
akan Nilsson.
Chalmers University of Technology,
Gothenburg, Sweden
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
December 10, 2013
1 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Overview
Ulf Nilsson
I
Introduction
I
Theory
I
Implementation
I
Modifications
I
Results
adjointShapeOptimizationFoam
December 10, 2013
2 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Introduction
Optimization in CFD
I
I
Increased importance
Optimization algorithm
-
I
Simplex
Evolution Strategies
Gradient based
...
Demanding
- Number of variables
- Constraints
Ulf Nilsson
adjointShapeOptimizationFoam
December 10, 2013
3 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Introduction
Adjoint approach
Ulf Nilsson
I
Gradient based method
I
Tool to compute sensitivity w.r.t. design variables
I
Well behaved objective function
I
Large number of design variables
I
Constraints, number of functions
adjointShapeOptimizationFoam
December 10, 2013
4 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Optimization problem
Objective function
minimize J = J(v, p, α)
(1)
such that
(v · ∇)v + ∇p − ∇ · (2νD(v)) + αv
= 0,
(2)
∇·v
= 0.
(3)
I
Gerneral cost function
I
Incompressible, steady state RANS equation
I
Porosity α
I
Topology optimization
Ulf Nilsson
adjointShapeOptimizationFoam
December 10, 2013
5 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Optimization problem
Lagrangian relaxation
Z
minimize L := J +
(u, q)RdΩ.
(4)
Ω
I
Lagrange multipliers, (u, q)
I
A penalty to violate the constraints
I
Adjoint velocity, u, and pressure, q
I
No physical meaning
Ulf Nilsson
adjointShapeOptimizationFoam
December 10, 2013
6 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Optimization problem
Variations w.r.t. flow and design variables
δL = δα L + δv L + δp L,
(5)
δv L + δp L = 0.
I
I
I
I
(6)
Total variation.
Choose adjoint variables so that the variation with respect to flow
variables vanishes.
Equation of the adjoint variables and an expression for the sensitivity
of L with respect to the porosity of each cell, derivations out of the
scope of the project.
Approximations needed to arrive at the adjoint equations.
- Forzen turbulence.
- Ducted flow specialization.
I
Integration by parts to arrive at volume contribution and boundary
contribution.
Ulf Nilsson
adjointShapeOptimizationFoam
December 10, 2013
7 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Resulting equations
Sensitivity
∂L
= ui · vi Vi ,
∂αi
(7)
I
Expression for cell i
I
Vi the cell volume.
I
Equation of the adjoint variables and an expression for the sensitivity
of L with respect to the porosity of each cell
Ulf Nilsson
adjointShapeOptimizationFoam
December 10, 2013
8 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Resulting equations
Adjoint equations
−2D(u)v = −∇q + ∇ · (2νD(u)) − αu −
∇·u=
I
∂JΩ
,
∂v
∂JΩ
.
∂p
(8)
(9)
Contribution from volumetric part of the cost function, JΩ .
Ulf Nilsson
adjointShapeOptimizationFoam
December 10, 2013
9 / 40
Overview
Introduction
Theory
Implementation
Modifications
Results
Resulting equations
Adjoint boundary conditions
Wall and inlet boundary conditions
ut = 0,
(10)
∂JΩ
,
un = −
∂p
n · ∇q = 0.
(11)
(12)
Outlet boundary conditions
q = u · v + un vn + ν(n · ∇)un +
0 = vn ut + ν(n · ∇)ut +
∂JΓ
.
∂vt
∂JΓ
,
∂vn
(13)
(14)
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
10 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Resulting equations
Adjoint boundary conditions
v
p
Wall
No-slip
Zero gradient
Inlet
Prescribed value
Zero gradient
Outlet
Zero gradient
p=0
Table: Boundary conditions for the primal quantities v and p.
I
Only valid for specific primal boundary condition.
I
Contribution from both the volume and boundary part of the cost
function, JΩ and JΓ respectively.
I
Eq. 13 and 14 used to calculate the adjoint pressure and tangential
part of the adjoint velocity at the outlet.
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
11 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Resulting equations
Summary of the theory
I
Introducing a Lagrange relaxation, with multipliers u and q.
I
Choice of adjoint variables gives adjoint equations and an expression
for the sensititivity.
I
Lengthy derivations.
I
Approximations and requirements to ensure a valid method.
I
If the cost function satisfies Jω = 0 the dependence of J enters the
solver only through the boundary conditions.
I
The optimization problem reduces to solving two flow equations, the
primal and the adjoint.
I
When the sensitivity is calculated an algorithm is needed to update
the porosity.
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
12 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Steepest descent algorithm
Steepest descent
pk = −∇f (xk ),
(15)
αn+1 = αn − ui · vi Vi δ.
(16)
I
General search direction in the opposite direction of the gradient of
the objective function.
I
Considering a linear problem a minimum will be found for small
enough steps in the search direction.
I
Using the sensitivity, ∂L/∂α to update.
I
The step length δ. An optimization problem of its own.
I
Under relaxation and limits in the final implementation.
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
13 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
Solution procedure
Main-loop
Solve the
primal system
(v, p)
“Frozen tubulence”
Update the
porosity
Solve the
adjoint system
Steepest descent
(u, q)
Convergence
assessment
Converged
End time or
End
convergence condition
Figure 1: Block scheme of the solution procedure used in adjointShapeOptimizationFoam.
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
14 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
$FOAM SOLVERS/incompressible/
adjointShapeOptimizationFoam/
adjointShapeOptimizationFoam.C
Let us study the implementation together, trying to identify the underlying
theory we have covered so far. In the order it is implemented in the code.
I
Porosity update
- Implementation of the steepest descent using under relaxation.
- Correct?
I
The primal pressure-velocity SIMPLE corrector
- How is the sink term implemented?
I
The adjoint pressure-velocity SIMPLE corrector
- Is the current implementation done for objective functions with
volumteric contribution (JΩ 6= 0)?
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
15 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
Porosity update
94
95
96
a l p h a +=
mesh . f i e l d R e l a x a t i o n F a c t o r ( " alpha " )
∗ ( min ( max ( a l p h a + lambda ∗ ( Ua & U) , z e r o A l p h a ) , alphaMax ) − a l p h a ) ;
Listing 1: file adjointShapeOptimizationFoam.C
αn+1 = αn (1 − γ) + γmin (max ((αn + ui · vi Vi δ) .0) , αmax ) ,
(17)
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
16 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
Porosity update
However the there seems to be an incorrect sign in the steepest descent
implementation. The correct eq. and implementation should be
αn+1 = αn (1 − γ) + γmin (max ((αn − ui · vi Vi δ) , 0) , αmax ) .
94
95
96
(18)
a l p h a +=
mesh . f i e l d R e l a x a t i o n F a c t o r ( " alpha " )
∗ ( min ( max ( a l p h a − lambda ∗ ( Ua & U) , z e r o A l p h a ) , alphaMax ) − a l p h a ) ;
Listing 2: file adjointShapeOptimizationFoam.C
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
17 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
Primal SIMPLE loop
103
104
105
106
107
108
109
110
111
112
113
114
// Momentum p r e d i c t o r
tmp<f v V e c t o r M a t r i x > UEqn
(
fvm : : d i v ( p h i , U)
+ t u r b u l e n c e −>d i v D e v R e f f (U)
+ fvm : : Sp ( a l p h a , U)
);
UEqn ( ) . r e l a x ( ) ;
s o l v e ( UEqn ( ) == −f v c : : g r a d ( p ) ) ;
Listing 3: file adjointShapeOptimizationFoam.C
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
18 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
Primal SIMPLE loop
I
fvm::Sp(alpha, U) Implicit implementation of the extra
source term.
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
19 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
Adjoint SIMPLE loop
156 // A d j o i n t Momentum p r e d i c t o r
157
158
v o l V e c t o r F i e l d a d j o i n t T r a n s p o s e C o n v e c t i o n ( ( f v c : : g r a d ( Ua ) & U ) ) ;
159
// v o l V e c t o r F i e l d a d j o i n t T r a n s p o s e C o n v e c t i o n
160
// (
161
// f v c : : r e c o n s t r u c t
162
// (
163
//
mesh . magSf ( ) ∗ ( f v c : : sn Grad ( Ua ) & f v c : : i n t e r p o l a t e (U) )
164
// )
165
// ) ;
166
167
zeroCells ( adjointTransposeConvection , i n l e t C e l l s );
168
169
tmp<f v V e c t o r M a t r i x > UaEqn
170
(
171
fvm : : d i v (− p h i , Ua )
172
− adjointTransposeConvection
173
+ t u r b u l e n c e −>d i v D e v R e f f ( Ua )
174
+ fvm : : Sp ( a l p h a , Ua )
175
);
176
177
UaEqn ( ) . r e l a x ( ) ;
178
179
s o l v e ( UaEqn ( ) == −f v c : : g r a d ( pa ) ) ;
Listing 4: file adjointShapeOptimizationFoam.C
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
20 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Solver
Adjoint SIMPLE loop
I
Similar to the primal system.
I
No additional terms containing information about the cost
function, i.e. current implementation treats cost functions of
the type JΓ = 0
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
21 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Adjoint boundary conditions
The bad results when the sign in the steepest descent algorithm is
changed indicates additional problems with the current solver. Since the
actual solver seems to agree with the theory, what remains to be examined
is the boundary conditions.
First the boundary condition for a specific cost function need to be
calculated.
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
22 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Power dissipation
Z
J := −
1
(p + v 2 )v · n dΓ.
2
(19)
Γ
I
Only boundary contribution.
I
Energy loss, due to pressure drop and kinetic energy.
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
23 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Inlet boundary conditions
ut = 0,
un = −
(20)
∂JΩ
,
∂p
n · ∇q = 0.


ut




un




n · ∇q
=(
0,
=
0
vn
(21)
(22)
at wall,
at inlet,
(23)
= 0.
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
24 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Inlet boundary conditions
I
Use existing conditions
- Prescribed value
- No slip
- Zero-gradient
I
Could be useful to implement a inlet boundary condition for
the adjoint velocity, if the velocity profile of the primal velocity
at the inlet is not easily described.
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
25 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Outlet boundary conditions
q = u · v + un vn + ν(n · ∇)un +
0 = vn ut + ν(n · ∇)ut +
(
q
0
∂JΓ
,
∂vn
(24)
∂JΓ
.
∂vt
(25)
= u · v + un vn + ν(n · ∇)un − 21 v 2 − vn2 ,
= vn (ut − vt ) + ν(n · ∇)ut .
(26)
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
26 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Outlet boundary conditions
I
Use eq. 26 to prescribe a value to the adjoint pressure
I
Use eq. 26 to prescribe a value to the tangential component
of the adjoint velocity
ut =
ν
−∆
uneigh,t − vp,n vp,t
,
ν
vp,n + ∆
(27)
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
27 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Implementation of pressure condition
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
// ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ Member F u n c t i o n s
∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ //
v o i d Foam : : a d j o i n t O u t l e t P r e s s u r e F v P a t c h S c a l a r F i e l d : : u p d a t e C o e f f s ( )
{
i f ( updated ( ) )
{
return ;
}
const f v s P a t c h F i e l d <s c a l a r >& p h i p =
p a t c h ( ) . l o o k u p P a t c h F i e l d <s u r f a c e S c a l a r F i e l d , s c a l a r >("phi" ) ;
const f v s P a t c h F i e l d <s c a l a r >& p h i a p =
p a t c h ( ) . l o o k u p P a t c h F i e l d <s u r f a c e S c a l a r F i e l d , s c a l a r >("phia" ) ;
const f v P a t c h F i e l d <v e c t o r >& Up =
p a t c h ( ) . l o o k u p P a t c h F i e l d <v o l V e c t o r F i e l d , v e c t o r >("U" ) ;
const f v P a t c h F i e l d <v e c t o r >& Uap =
p a t c h ( ) . l o o k u p P a t c h F i e l d <v o l V e c t o r F i e l d , v e c t o r >("Ua" ) ;
o p e r a t o r ==(( p h i a p / p a t c h ( ) . magSf ( ) − 1 . 0 ) ∗ p h i p / p a t c h ( ) . magSf ( ) + ( Up & Uap ) ) ;
fixedValueFvPatchScalarField : : updateCoeffs ( ) ;
}
Listing 5: file adjointOutletPressureFvPatchScalarField.C
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
28 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Implementation of velocity condition
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
// ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ Member F u n c t i o n s
∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ ∗ //
// Update t h e c o e f f i c i e n t s a s s o c i a t e d w i t h t h e p a t c h f i e l d
v o i d Foam : : a d j o i n t O u t l e t V e l o c i t y F v P a t c h V e c t o r F i e l d : : u p d a t e C o e f f s ( )
{
i f ( updated ( ) )
{
return ;
}
const f v s P a t c h F i e l d <s c a l a r >& p h i a p =
p a t c h ( ) . l o o k u p P a t c h F i e l d <s u r f a c e S c a l a r F i e l d , s c a l a r >("phia" ) ;
const f v P a t c h F i e l d <v e c t o r >& Up =
p a t c h ( ) . l o o k u p P a t c h F i e l d <v o l V e c t o r F i e l d , v e c t o r >("U" ) ;
s c a l a r F i e l d Un ( mag ( p a t c h ( ) . n f ( ) & Up ) ) ;
v e c t o r F i e l d UtHat ( ( Up − p a t c h ( ) . n f ( ) ∗ Un ) / ( Un + SMALL ) ) ;
v e c t o r F i e l d Uan ( p a t c h ( ) . n f ( ) ∗ ( p a t c h ( ) . n f ( ) & p a t c h I n t e r n a l F i e l d ( ) ) ) ;
v e c t o r F i e l d : : o p e r a t o r =( p h i a p ∗ p a t c h ( ) . S f ( ) / s q r ( p a t c h ( ) . magSf ( ) ) + UtHat ) ;
// v e c t o r F i e l d : : o p e r a t o r =(Uan + UtHat ) ;
fixedValueFvPatchVectorField : : updateCoeffs ( ) ;
}
Listing 6: file adjointOutletVelocityFvPatchVectorField.C
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
29 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Analyzing the implementation
Listing 5 gives an implementation of the pressure condition according to
q = (un − 1)vn + u · v.
(28)
1
q = u · v + un vn + ν(n · ∇)un − v 2 − vn2 ,
2
(29)
compare to
Listing 6 gives an implementation of the velocity condition according to
up,t =
vp − vp,n
.
up,n + SMALL
(30)
compare to
ut =
ν
−∆
uneigh,t − vp,n vp,t
,
ν
vp,n + ∆
(31)
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
30 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Boundary conditions of the adjoint variables
Analyzing the theory
I
Not equal to the theory.
I
Another cost function, total pressure loss.
I
Approximations or derivations done makes the cost function
hard to identify.
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
31 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Analyzing the theory
I
Changing the sign of the steepest descent implementation.
I
Implementing the cost functions according to derived
equations.
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
32 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Adjoint pressure boundary condition
I
I
Similar to the existing implementation
A few additional terms
- Velocity in the neighbouring node
- Effective viscosity
const s c a l a r F i e l d & d e l t a i n v = p a t c h ( ) . d e l t a C o e f f s ( ) ; // d i s t a n c e ˆ( −1) between p a t c h and n e i g h b o u r i n g node .
// c r e a t e t h e o b j e c t n e e d e d t o g e t t h e v i s c o u s i t y :
const i n c o m p r e s s i b l e : : RASModel& r a s M o d e l =
db ( ) . l o o k u p O b j e c t <i n c o m p r e s s i b l e : : RASModel>(" RASProperties " ) ;
// n u e f f from t h e t u r u l e n c e model , r a s M o d e l :
s c a l a r F i e l d n u e f f = rasModel . nuEff ( ) ( ) . boundaryField ( ) [ patch ( ) . index ( ) ] ;
// N e i g h b o u r i n g node ’ s v e l o c i t y ( n o r m a l component ) :
s c a l a r F i e l d U a n e i g h n = ( Uap . p a t c h I n t e r n a l F i e l d ( ) & p a t c h ( ) . n f ( ) ) ;
Listing 7: file myAdjointOutletVelocityFvPatchVectorField.C
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
33 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Adjoint velocity boundary condition
I
Vector fields instead of scalar fields
// p a t c h I n t e r n a l F i e l d t o g e t t h e a d j o i n t v e l o c i t y o f n e i g h b o u r i n g node .
v e c t o r F i e l d U a n e i g h = Uap . p a t c h I n t e r n a l F i e l d ( ) ;
// I t s t a n g e n t i a l and n o r m a l components
v e c t o r F i e l d Uaneigh n = ( Uaneigh & normal )∗ normal ;
v e c t o r F i e l d Uaneigh t = Uaneigh − Uaneigh n ;
Listing 8: file myAdjointOutletVelocityFvPatchVectorField.C
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
34 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Case settings
The boundary conditions of the primal and adjoint variables are set
according to the table below. More detailed information and m4 script of
the “box” example case can be found in the provided files.
u
q
Wall
Prescribed value, 0
Zero gradient
Inlet
Prescribed value, un = vn
Zero gradient
Outlet
adjointOutletVelocityPower
adjointOutletPressurePower
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
35 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Power dissipation
pitzDaily porosity field
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
36 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Power dissipation
pitzDaily value of the objective function
0.006
”valueObjFunc.xy”
0.005
J
0.004
0.003
0.002
0.001
0
0
500
1000
1500
2000
2500
3000
Iteration step
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
37 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Power dissipation
pitzDaily value of the objective function, without
porosity update
0.006
”valueObjFunc.xy”
0.005
J
0.004
0.003
0.002
0.001
0
0
500
1000
1500
2000
2500
3000
Iteration step
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
38 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Power dissipation
Pipe bend example
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
39 /
40
Overview
Introduction
Theory
Implementation
Modifications
Results
Power dissipation
Thank you for listening!
Questions
December 10, 2013
Ulf Nilsson
adjointShapeOptimizationFoam
40 /
40