Table Of ContentOptimal Reconstruction of Material Properties in
Complex Multiphysics Phenomena
Vladislav Bukshtynov and Bartosz Protas
∗ †
3
1
0
2 ∗School of Computational Science and Engineering, McMaster University
n 1280 Main Street West, Hamilton, Ontario, CANADA L8S 4K1
a
e-mail: [email protected]
J
4
2 †Department of Mathematics and Statistics, McMaster University
1280 Main Street West, Hamilton, Ontario, CANADA L8S 4K1
]
n
e-mail: [email protected]
y
d
- Abstract
u
fl We develop an optimization–based approach to the problem of reconstructing
.
s temperature–dependentmaterialpropertiesincomplexthermo–fluidsystemsdescribed
c
by the equations for the conservation of mass, momentum and energy. Our goal is to
i
s
estimate the temperature dependence of the viscosity coefficient in the momentum
y
h equation based on some noisy temperature measurements, where the temperature is
p
governed by a separate energy equation. We show that an elegant and computation-
[
ally efficient solution of this inverse problem is obtained by formulating it as a PDE–
1
constrained optimization problem which can be solved with a gradient–based descent
v
2 method. A key element of the proposed approach, the cost functional gradients are
2 characterized by mathematical structure quite different than in typical problems of
9
PDE–constrained optimization and are expressed in terms of integrals defined over the
5
. level sets of the temperature field. Advanced techniques of integration on manifolds
1
are required to evaluate numerically such gradients, and we systematically compare
0
3 three different methods. As a model system we consider a two–dimensional unsteady
1
flow in a lid–driven cavity with heat transfer, and present a number of computational
:
v tests to validate our approach and illustrate its performance.
i
X
r Keywords: parameter estimation, material properties, optimization, adjoint anal-
a
ysis, integration on level sets.
1 Introduction
In this work we propose and validate a computational approach to the reconstruction of
material properties in complex multiphysics phenomena based on incomplete and possibly
1
noisy measurements. The material properties we are interested in here are the transport co-
efficients characterizing diffusion processes such as the viscosity or the thermal conductivity,
and we focus on problems in which these coefficients depend on the state variables in the
system. By the “multiphysics” aspect we mean situations in which the material property
used in one conservation equation is a function of a state variable governed by a different
conservation equation, e.g., reconstruction of the temperature dependence of the viscosity
coefficient used in the momentum equation, where the temperature is governed by a sepa-
rate energy equation, which is the specific model problem investigated in this study. This
research is motivated by questions arising in the computational analysis and optimization of
advanced welding processes which involves modelling complex alloys in the liquid phase at
high temperatures [1].
Inverse problems of parameter estimation for partial differential equations (PDEs) have
received significant attention in the literature, both regarding theoretical [2] and practical
aspects [3]. However, most of the attention focused on problems in which the material prop-
erties are functions of the space variable (i.e., the independent variable in the problem).
Such problems are, at least in principle, relatively well understood and represent the foun-
dation of, for example, numerous imaging techniques in medicine [4] and earth sciences [5].
The problem considered here is in fact different, in that the material properties are sought
as functions of the state (dependent) variables in the system which gives rise to a number
of computational challenges absent in the “classical” parameter estimation problem. Other
than the seminal work of Chavent and Lemonnier [6], our earlier study [7] concerning a
simplified model problem and a few investigations of fully discrete formulations surveyed
in [7], there does not seem to be much literature concerning computational methods for
this type of parameter estimation problems. One way to solve such inverse problems is to
formulate them as optimization problems and this is the approach we will follow focusing on
the “optimize–then–discretize” paradigm in which the optimality conditions are formulated
at the continuous (PDE) level and only then discretized. The goal of this investigation is to
extend the approach formulated in [7] for a simple model to a realistic multiphysics problem
involving time–dependent fluid flow in a two–dimensional (2D) domain. As will be shown
below, a number of computational difficulties will need to be overcome in order to achieve
this goal.
As a key contribution of this work, we address a number of computational challenges
related to accurate and efficient evaluation of cost functional gradients which are critical to
the implementation of the proposed approach. More specifically, these gradients are given in
terms of integrals of expressions involving state and adjoint variables defined on a grid over
contours given by the level sets of the temperature field. A number of techniques have been
proposed for the numerical evaluation of integrals defined over manifolds defined by level–set
functions. Some of them rely on regularized Dirac delta and Heaviside functions [8,9], or
discretization of the Dirac delta function [10–12]. Similar approaches, based on approxima-
tions of the Dirac delta functions obtained using the level–set function and its gradient, were
developed by Towers [13,14]. The family of geometric approaches proposed by Min and Gi-
bou in [15,16] relies on a decomposition of the entire domain into simplices. We emphasize
2
that the problem discussed in this work is in fact more complicated, as the computation
of our cost functional gradients requires evaluation of the corresponding integrals for the
level–set values spanning the entire state space of interest, hence there are also additional
issues related to the discretization of the state space which were outside the scope of refer-
ences [8–16]. Thus, in order to address these questions and assess the different trade–offs in
the choice of the numerical parameters we will compare the computational performance of
three different methods for the evaluation of cost functional gradients.
The structure of this paper is as follows: in the next Section we formulate our model
problem, in the following Section we cast the problem of parameter estimation as an opti-
mization problem, an adjoint–based gradient–descent algorithm is formulated in Section 4,
in Section 5 we outline some regularization strategies needed in the presence of measurement
noise, whereas in Section 6 we analyze three different numerical approaches to the evaluation
of the cost functional gradients; extensive computational results are presented in Section 7
with discussion and conclusions deferred to Section 8.
2 Model Problem
Let Ω Rd, d = 2,3, be the spatial domain on which our model problem is formulated.
⊂
To fix attention, but without loss of generality, in the present investigation we focus on
the problem of a reconstruction of the temperature dependence µ = µ(T) of the viscosity
coefficient µ : R R+ in the momentum equation (Navier-Stokes equation), where the
→
temperature T is governed by a separate energy equation
(cid:2) (cid:3)
∂ u+(u ∇)u+∇p ∇ µ(T)[∇u+(∇u)T] =0 in Ω, (1a)
t
· − ·
∇ u =0 in Ω, (1b)
·
∂ T +(u ∇)T ∇ [k∇T] =0 in Ω, (1c)
t
· − ·
subject to appropriate Dirichlet (or Neumann) boundary and initial conditions
u = u on ∂Ω, (2a)
B
T = T on ∂Ω, (2b)
B
u( ,0) = u , T( ,0) = T in Ω. (2c)
0 0
· ·
The specific inverse problem we address in this investigation is formulated as follows. Given
a set of time–dependent “measurements” T˜(t) M of the state variable (temperature) T
{ i }i=1
at a number of points x M in the domain Ω (or along the boundary ∂Ω) and obtained
{ i}i=1
within the time window t [0,t ], we seek to reconstruct the constitutive relation µ = µ(T)
f
∈
such that solutions of problem (1)-(2) obtained with this reconstructed function will fit best
the available measurements.
In regard to reconstruction of constitutive relations in general, it is important that such
relations be consistent with the second principle of thermodynamics [17]. There exist two
mathematical formalisms, one due to Coleman and Noll [18] and another one due to Liu [19],
3
developed to ensure in a very general setting that a given form of the constitutive relation
does not violate the second principle of thermodynamics. In continuous thermodynamical
andmechanicalsystemsthisprincipleisexpressedintermsoftheClausius–Duheminequality
[20] which in the case of our present model problem (1)–(2) reduces to the statement that
µ(T) > 0 for all possible values of T.
In our discussion below we will also need definitions of the following intervals, cf. Figure
1:
[T ,T ] (cid:44) [min T(x),max T(x)]whichrepresentsthetemperaturerangespanned
• α β x∈Ω x∈Ω
by the solution of problem (1); thus, following [21], we will refer to the interval
(cid:44) [T ,T ] as the identifiability interval,
α β
I
L(cid:32) (cid:44) [T ,T ], where T T and T T ; this will be the temperature interval on
a b a α b β
• ≤ ≥
which we will seek to obtain a reconstruction of the material property; we note that in
general the interval L(cid:32) will be larger than the identifiability interval, i.e., L(cid:32) , and
I ⊆
(cid:44) [min min T˜(t),max max T˜(t)] which defines the temper-
1≤i≤M 0<t≤t i 1≤i≤M 0<t≤t i
• M f f
ature range spanned by the measurements T˜ M ; this interval is always contained
{ i}i=1
the identifiability interval , i.e., ; we will refer to the interval as the
I M ⊆ I M
measurement span.
3 Parameter Estimation as an Optimization Problem
It is assumed that the constitutive relations µ(T) are differentiable functions of the state
variable (temperature) and belong to the following set
= µ(T) piecewise C1 on L(cid:32) ; 0 < m µ(T) M < , T L(cid:32) , (3)
µ µ µ
S { ≤ ≤ ∞ ∀ ∈ }
where m ,M R+. We will also assume that the set consisting of constitutive relations
µ µ µ
∈ S
µ(T) defined on L(cid:32) is embedded in a Hilbert (function) space to be specified below. Solving
X
ourparameterestimationproblemisthereforeequivalenttofindingasolutiontotheoperator
equation
(µ) = T, (4)
F
where : (L ([0,t ]))M is the map from the constitutive relations to the measure-
µ 2 f
F S →
ments. An approach commonly used to solve such problems consists in reformulating them
as least–squares minimization problems which in the present case can be done by defining
the cost functional : R as
J X →
1 (cid:90) tf (cid:88)M (cid:104) (cid:105)2
(µ) (cid:44) T(τ,x ;µ) T˜(τ) dτ, (5)
i i
J 2 −
0
i=1
4
Figure 1: Schematic showing (left) the solution T(t ,x) at some fixed time t and (right)
0 0
the corresponding constitutive relation µ(T) defined over their respective domains, i.e., Ω =
( 1,1) and the identifiability region . The thick dotted line represents an extension of the
− I
constitutive relation µ(T) from to the interval L(cid:32) . In the Figure on the right the horizontal
I
axis is to be interpreted as the ordinate.
5
where the dependence of the temperature field T( ;µ) on the form of the constitutive relation
·
µ = µ(T) is given by governing system (1)–(2). The optimal reconstruction µˆ is obtained as
a minimizer of cost functional (5), i.e.,
µˆ = argmin (µ). (6)
J
µ∈Sµ
We recall that the constitutive property is required to satisfy the positivity condition µ(T) >
0 for all T L(cid:32) , cf. (3). Therefore, the optimal reconstruction µˆ should in fact be obtained
∈
as an inequality–constrained minimizer of cost functional (5), i.e.,
µˆ = argmin (µ). (7)
J
µ∈X,
L(cid:32)
µ(T)>0, T∈
We add that in problems involving constitutive relations depending on several state vari-
ables the inequality constraint µ(T) > 0 will be replaced with a more general form of
the Clausius–Duhem inequality [20]. Different computational approaches for converting
inequality–constrained optimization problems to unconstrained formulations are surveyed
in [22,23]. Here we follow a straightforward approach based on the so-called “slack” vari-
able [24]. We define a new function θ(T) : R R such that
→
µ(T) = θ2(T)+m , (8)
µ
where m is a lower bound for µ(T), cf. (3). This change of variables allows us to transform
µ
the inequality–constrained optimization problem (7) to a new unconstrained one
ˆ
θ = argmin (θ), (9)
J
θ∈X
where the constraint µ(T) > 0 is satisfied automatically when minimization is performed
with respect to the new variable θ(T). In view of (3), we note that the new optimization
variable θ belongs to the following set
(cid:112)
= θ(T) piecewise C1 on L(cid:32) ; θ(T) < M m , T L(cid:32) . (10)
θ µ µ
S { | | − ∀ ∈ }
The governing PDE system (1) can thus be rewritten in the form
(cid:2) (cid:3)
∂ u+(u ∇)u+∇p ∇ (θ2(T)+m )[∇u+(∇u)T] =0 in Ω, (11a)
t µ
· − ·
∇ u =0 in Ω, (11b)
·
∂ T +(u ∇)T ∇ [k∇T] =0 in Ω, (11c)
t
· − ·
subject to Dirichlet boundary and initial conditions (2). The new problem (9) requires
redefining cost functional (5) in terms of the new variable
1 (cid:90) tf (cid:88)M (cid:104) (cid:105)2
˜
(θ) = T(τ,x ;θ) T (τ) dτ. (12)
i i
J 2 −
0
i=1
6
Problem (9) is characterized by the first–order optimality condition which requires the
Gaˆteaux differential of cost functional (12), defined as (cid:48)(θ;θ(cid:48)) (cid:44) lim (cid:15)−1[ (θ + (cid:15)θ(cid:48))
(cid:15)→0
J J −
(θ)], to vanish for all perturbations θ(cid:48) [25], i.e.,
J ∈ X
(cid:48)(θˆ;θ(cid:48)) = 0. (13)
θ(cid:48)∈X
∀ J
ˆ
The (local) optimizer θ can be computed with the following gradient descent algorithm as
θˆ= lim θ(n), where
n→∞
(cid:40)
θ(n+1) = θ(n) τ(n)∇ (θ(n)), n = 1,...,
θ
− J (14)
θ(1) = θ ,
0
in which ∇ (θ) represents the gradient of cost functional (θ) with respect to the control
θ
J J
variableθ (wewilladopttheconventionthatasubscriptontheoperator∇willbeusedwhen
differentiation is performed with respect to variables other than x), τ(n) is the length of the
step along the descent direction at the n–th iteration, whereas θ = √µ m is the initial
0 0 µ
−
guess taken, for instance, corresponding to a constant µ , or some other appropriate initial
0
approximation. For the sake of clarity, formulation (14) represents the steepest–descent
algorithm, however, in practice one typically uses more advanced minimization techniques,
such as the conjugate gradient method, or one of the quasi–Newton techniques [26]. We note
that, since minimization problem (9) is in general nonconvex, condition (13) characterizes
only a local, rather than global, minimizer.
The key ingredient of minimization algorithm (14) is computation of the cost functional
gradient ∇ (θ). We emphasize that, since θ = θ(T) is a continuous variable, the gradient
θ
J
∇ (θ) represents in fact the infinite–dimensional sensitivity of (θ) to perturbations of
θ
J J
θ(T). This gradient can be determined based on suitably defined adjoint variables (Lagrange
multipliers) obtained from the solution of the corresponding adjoint system. Since this
derivation differs in a number of imported technical details from analogous derivations in
“standard” PDE–constrained optimization problems, it will be reviewed in Section 4. The
expression for the gradient is then validated for consistency in Section 7.4.
4 Cost Functional Gradients via Adjoint–based Anal-
ysis
Since the new variable θ(T) belongs to set , cf. (10), we will seek to reconstruct θ(T) as
θ
S
elements of the Sobolev space H1(L(cid:32) ), so that the gradient ∇ will need to be obtained
θ
J
with respect to the corresponding inner product. However, in order to make the derivation
procedure easier to follow, we will first obtain an expression for the gradient in the space
L (L(cid:32) ), and only then will obtain the Sobolev gradients which will be eventually used in the
2
solution of optimization problem (9). In all these steps our transformations will be formal.
We begin by computing the directional (Gaˆteaux) differential of cost functional (12) which
7
yields
(cid:90) t M
f (cid:88)
(cid:48)(θ;θ(cid:48)) = [T(τ,x ;θ) T˜(τ)]T(cid:48)(τ,x ;θ,θ(cid:48))dτ, (15)
i i i
J −
0
i=1
where the perturbation variable T(cid:48)(θ,θ(cid:48)) satisfies the perturbation system obtained from
(1)–(2). Next, we invoke the Riesz representation theorem [27] for the directional differential
(cid:48)(θ; ) which yields
J · (cid:68) (cid:69)
(cid:48)(θ;θ(cid:48)) = ∇ ,θ(cid:48) , (16)
θ
J J X
where , represents an inner product in the Hilbert space (we will first set = L (L(cid:32) )
X 2
(cid:104)· ·(cid:105) X X
and afterwards change this to = H1(L(cid:32) )). We note that the expression on the right–hand
X
side (RHS) in (15) is not consistent with Riesz representation (16), since, as will be shown
below, the perturbation variable θ(cid:48) is hidden in the system defining T(cid:48)(θ,θ(cid:48)). However, this
expression can be transformed to Riesz form (16) with the help of a suitably–defined adjoint
variable, a result which is stated in Theorem 4.1 below. The main aspect in which this
derivation differs from the standard adjoint analysis [28] is that the inner product in Riesz
identity (16) is defined using the state variable (temperature) as the integration variable,
whereas the variational formulation of (1)–(2) is defined using integration with respect to
the independent variables (x and t).
Theorem 4.1. Let Ω be a sufficiently regular open bounded domain and θ(cid:48) = H1(L(cid:32) ).
∈ X
We assume that the solutions u and T of system (1)–(2) are sufficiently smooth. Then, the
Riesz representation of directional differential (15) has the form
(cid:90) ∞ (cid:90) (cid:20)(cid:90) t (cid:21)
f
(cid:48)(θ;θ(cid:48)) = 2 δ(T(x) s)θ(s) [∇u+(∇u)T] : ∇u∗dτ θ(cid:48)(s)dxds, (17)
J − −
−∞ Ω 0
where δ( ) denotes Dirac delta function and the adjoint state u∗,T∗ is defined as the
· { }
solution of the system
∂ u∗ (u ∇)u∗ ∇ σ∗ +u∗ (∇u)T +T∗∇T = 0 in Ω,
t
− − · − · ·
(18a)
∇ u∗ = 0 in Ω,
·
(18b)
dθ
∂ T∗ (u ∇)T∗ ∇ [k∇T∗]+2θ(T) (T)[∇u+(∇u)T] : ∇∗u
t
− − · − · dT
M
(cid:88)
˜
= [T(x ;θ) T ]δ(x x ) in Ω,
i i i
− −
i=1
(18c)
where σ∗ (cid:44) p∗ +(θ2(T)+m )(cid:2)∇u∗ +(∇u∗)T(cid:3), with the following boundary and terminal
µ
− I
8
conditions
u∗ = 0 on ∂Ω,
T∗ = 0 on ∂Ω, (19)
u∗( ;t ) = 0, T∗( ;t ) = 0 in Ω.
f f
· ·
Proof. We will denote the stress tensor σ (cid:44) p +(θ2(T)+m )(cid:2)∇u+(∇u)T](cid:3) and rewrite
µ
− I
the governing system (11) as
∂ u+(u ∇)u ∇ σ = 0 in Ω,
t
· − ·
∇ u = 0 in Ω, (20)
·
∂ T +(u ∇)T ∇ [k∇T] = 0 in Ω
t
· − ·
with the boundary and initial conditions (2). Perturbing the state variables u, p and T,
which are functions of time and space, we get
u = u +(cid:15)u(cid:48) +Ø((cid:15)2),
0
p = p +(cid:15)p(cid:48) +Ø((cid:15)2), (21)
0
T = T +(cid:15)T(cid:48) +Ø((cid:15)2),
0
so that the corresponding expansion of the constitutive relation θ(T) will have the following
form
dθ
θ(T) = θ (T)+(cid:15)θ(cid:48)(T)+Ø((cid:15)2) = θ (T )+(cid:15) (T )T(cid:48)(θ ;θ(cid:48))+(cid:15)θ(cid:48)(T )+Ø((cid:15)2), (22)
0 0 0 0 0 0
dT
where the subscript “0” is used to denote the unperturbed (reference) material property and
the state variable, whereas the prime denotes the corresponding perturbations. We also have
dθ
θ2(T) = θ2(T )+2(cid:15)θ (T ) (T )T(cid:48)(θ ;θ(cid:48))+2(cid:15)θ (T )θ(cid:48)(T )+Ø((cid:15)2). (23)
0 0 0 0 dT 0 0 0 0 0
Substituting (21) and (23) into (20), collecting terms corresponding to (cid:15) in different powers
and denoting
σˆ (cid:44) p(cid:48) +(θ2(T)+m )(cid:2)∇u(cid:48) +(∇u(cid:48))T(cid:3),
µ
− I
(cid:18) (cid:19)
dθ
σ˜ (cid:44) 2θ(T) (T)T(cid:48)(θ ;θ(cid:48))+2θ(T)θ(cid:48)(T) (cid:2)∇u+(∇u)T(cid:3),
0
dT
we obtain the perturbation (sensitivity) system corresponding to (11)
∂ u(cid:48) +(u(cid:48) ∇)u+(u ∇)u(cid:48) ∇ (σˆ +σ˜) =0 in Ω, (24a)
t
· · − ·
∇ u(cid:48) =0 in Ω, (24b)
·
∂ T(cid:48) +(u(cid:48) ∇)T +(u ∇)T(cid:48) ∇ [k∇T(cid:48)] =0 in Ω (24c)
t
· · − ·
9
with the following boundary and initial conditions
u(cid:48) = 0 on ∂Ω, (25a)
T(cid:48) = 0 on ∂Ω, (25b)
u(cid:48)( ,0) = 0, T(cid:48)( ,0) = 0 in Ω. (25c)
· ·
Then, integrating equation (24a) against u∗, equation (24b) against p∗, and equation (24c)
against T∗ over the space domain Ω and time [0,t ], integrating by parts and factorizing u(cid:48),
f
T(cid:48) and p(cid:48), we arrive at the following relation
(cid:90) t (cid:90)
f (cid:2) (cid:3)
∂ u∗ +u∗ (∇u)T (u ∇)u∗ ∇ σ∗ +T∗∇T u(cid:48)dxdτ
t
− · − · − · ·
0 Ω
(cid:90) t (cid:90)
f
(∇ u∗)p(cid:48)dxdτ
− ·
0 Ω
(cid:90) tf (cid:90) (cid:20) dθ (cid:21)
+ ∂ T∗ (u ∇)T∗ ∇ (k∇T∗)+2θ(T) (T)[∇u+(∇u)T] : ∇u∗ T(cid:48)dxdτ
t
− − · − · dT
0 Ω
(cid:90) t (cid:90)
f
+ 2θ(T)θ(cid:48)(T)[∇u+(∇u)T] : ∇u∗dxdτ = 0.
0 Ω
(26)
We now require that the adjoint variables u∗, p∗ and T∗ satisfy system (18)–(19). We also
note that owing to the judicious choice of the RHS term in (18c), the last term in relation
(26) is in fact equal to the directional differential (cid:48)(θ;θ(cid:48)), so that we have
J
(cid:90) t (cid:90)
f
(cid:48)(θ;θ(cid:48)) = 2 θ(T(x,τ))θ(cid:48)(T(x,τ))[∇u(x,τ)+(∇u(x,τ))T] : ∇u∗(x,τ)dxdτ,
J −
0 Ω
(27)
where, for emphasis, we indicated the integration variables as arguments of the state and
adjoint variables. We note that this expression is still not in Riesz form (16), where inte-
gration must be performed with respect to the state variable (temperature T). Thus, we
proceed to express for any given function f(T) its pointwise evaluation at T(x) through the
following integral transform. It is defined using a “change–of–variable” operator, denoted Π,
such that for given functions f : R R and T : Ω R, we have
→ →
(cid:90) +∞
f(T(x)) = δ(T(x) s)f(s)ds (cid:44) (Πf)(x). (28)
−
−∞
Using this transform to express f(T(x)) = θ(T(x,τ))θ(cid:48)(T(x,τ)) in (27) and changing the
order of integration (Fubini’s Theorem), we obtain expression (17) which is the required
Riesz representation (16) of directional differential (15).
We remark that we were able to prove an analogous result using a simpler approach based
on the Kirchhoff transform in [7], where both the constitutive relation and the state variable
were governed by the same equation (i.e., the problem was not of the “multiphysics” type).
10
