# Improve Counterexample Quality

Li Huang  Chair of Software Engineering
Schaffhausen Institute of Technology
Schaffhausen, Switzerland
<EMAIL_ADDRESS>Bertrand Meyer  Chair of Software Engineering
Schaffhausen Institute of Technology
Schaffhausen, Switzerland
<EMAIL_ADDRESS>Manuel Oriol  Chair of Quantum Software Engineering
Schaffhausen Institute of Technology
Schaffhausen, Switzerland
<EMAIL_ADDRESS>
# Improve Counterexample Quality for Program Verification

Li Huang  Chair of Software Engineering
Schaffhausen Institute of Technology
Schaffhausen, Switzerland
<EMAIL_ADDRESS>Bertrand Meyer  Chair of Software Engineering
Schaffhausen Institute of Technology
Schaffhausen, Switzerland
<EMAIL_ADDRESS>Manuel Oriol  Chair of Quantum Software Engineering
Schaffhausen Institute of Technology
Schaffhausen, Switzerland
<EMAIL_ADDRESS>
# Improve Counterexample Quality for Eiffel Program Verifier

Li Huang  Chair of Software Engineering
Schaffhausen Institute of Technology
Schaffhausen, Switzerland
<EMAIL_ADDRESS>Bertrand Meyer  Chair of Software Engineering
Schaffhausen Institute of Technology
Schaffhausen, Switzerland
<EMAIL_ADDRESS>Manuel Oriol  Chair of Quantum Software Engineering
Schaffhausen Institute of Technology
Schaffhausen, Switzerland
<EMAIL_ADDRESS>
# Improving Counterexample Quality from Failed Program Verification

Li Huang  Chair of Software Engineering
Schaffhausen Institute of Technology
Schaffhausen, Switzerland
<EMAIL_ADDRESS>Bertrand Meyer  Chair of Software Engineering
Schaffhausen Institute of Technology
Schaffhausen, Switzerland
<EMAIL_ADDRESS>Manuel Oriol  Chair of Quantum Software Engineering
Schaffhausen Institute of Technology
Schaffhausen, Switzerland
<EMAIL_ADDRESS>
###### Abstract

In software verification, a successful automated program proof is the ultimate
triumph. The road to such success is, however, paved with many failed proof
attempts. The message produced by the prover when a proof fails is often
obscure, making it very hard to know how to proceed further. The work reported
here attempts to help in such cases by providing immediately understandable
counterexamples. To this end, it introduces an approach called Counterexample
Extraction and Minimization (CEAM). When a proof fails, CEAM turns the
counterexample model generated by the prover into a a clearly understandable
version; it can in addition simplify the counterexamples further by minimizing
the integer values they contain. We have implemented the CEAM approach as an
extension to the AutoProof verifier and demonstrate its application to a
collection of examples.

###### Index Terms:

Program Verification, Counterexample, AutoProof, Boogie, SMT

##  I Introduction

Deductive program verification performs a rigorous analysis of the correctness
of programs with respect to their functional behavior, usually specified
formally by _contracts_ (such as pre- and postconditions, can class and loop
invariants). The approached has progressed in recent years thanks to the
development of powerful proof engines. In practice, however, verifying
industrial applications remains difficult. One of the obstacles is the lack of
intuitive feedback to understand the reasons why a verification attempt
failed. Although in many cases the underlying prover can provide a
counterexample containing some usable diagnostic information for debugging,
such a counterexample usually contains hundreds of difficult-to-interpret
lines. Another obstacle to usability is that integer values generated by the
prover for the counterexamples are often very large, and hence do not provide
programmers with an easy intuitive understanding of what is wrong.

This article presents the Counterexample Extraction and Minimization (CEAM)
approach for improving the quality of counterexamples produced when a proof
fails, and making them usable for identifying and correcting the underlying
bugs. We have implemented CEAM as an extension of the AutoProof environment
[1, 23], a static verification platform for contract-equipped Eiffel [15]
programs based on Hoare-style proofs. AutoProof relies on the Boogie proof
system [2, 11] and takes advantage of Boogie’s underlying SMT (Satisfiability
Modulo Theories) solver, by default Z3 [7].

When a proof fails, CEAM exploits the counterexample models (hereinafter
referred to simply as models) generated by the SMT solver and generates simple
counterexamples in a format more intuitive to programmers. CEAM also provides
a minimization mechanism allowing programmers to get simplified
counterexamples with integer variables reduced to their minimal possible
values. The current version of CEAM supports primitive types (integer,
boolean), user-defined types (classes) as well as some commonly used container
types (such as arrays and sequences).

Section II illustrates an example of using the CEAM approach. Section III
introduces the technologies used in our verification framework. Section IV
describes the details of the implementation of the CEAM. We evaluate the
applicability of the CEAM through a series of examples in Section V. After a
review of related work in Section VI, Section VII concludes the paper with our
ongoing work.

##  II An Example Session

Before exploring the principles and technologies of the CEAM approach, we look
at its practical use on a representative example (Fig. 1). The intent of the
max function in class MAX is to return into Result the maximum element of an
integer array a of size a.count. The two postconditions in lines 22 and 23
(labeled by _is_max_ and _in_array_) specify this intent: every element of the
array should be less than or equal to Result; at least one element should be
equal to Result.

When we try to verify the max function using AutoProof, verification fails and
AutoProof returns an error message “Postcondition _is_max_ may be violated”
(the first row in Fig. 2). Such a generic message tells us that the prover
cannot establish the postcondition, but does not enable us to find out why. In
this case, programmers can look at the model generated by the Z3 solver to
understand the cause of the failure. Deciphering the model is a cumbersome
task: the model spans hundreds of lines and is expressed in a cryptic
formalism.

In contrast, AutoProof extended with CEAM automatically generates a much
simpler counterexample from the model. As displayed in the second row of Fig.
2, the counterexample contains the initial values (on entry of max) of the
array’s size and of some of its elements. Seeing these concrete values, rather
than just the prover’s general failure message, helps the programmer
conjecture possible reasons for the failure. The values in the counterexample
are large, however, too large to give the programmer a direct intuition of the
problem at a human scale.

⬇

1class MAX feature

2max (a: ARRAY $[$INTEGER$]$): INTEGER

3 require a.count $>$ $0$

4 local i: INTEGER

5 do

6 from

7 Result $:=$ a $[1]$; i $:=$ $2$

8 invariant

9 $2$ $\ \\!\leq$ i and i $\ \\!\leq$ a.count $+$ $1$

10 $\forall$ j: $1$ $|..|$ $($i - $1)$ | a.sequence $[$j$]$ $<=$ Result

11 $\exists$ j: $1$ $|..|$ $($i - $1)$ | a.sequence $[$j$]$ $=$ Result

12 until

13 i $\geq$ a.count

14 loop

15 if a $[$i$]$ $>\>$ Result then

16 Result $:=$ a $[$i$]$

17 end

18 i $:=$ i $+$ $1$

19 variant a.count - i

20 end

21 ensure

22 is_max: $\forall$ j: $1$ $|..|$ a.count | a.sequence $[$j$]$ $<=$ Result

23 in_array: $\exists$ j: $1$ $|..|$ a.count | a.sequence $[$j$]$ $=$ Result

24 end

25end

Figure 1: MAX is a class that finds the maximum element of an integer array; a
fault (the exit condition at line 13 is incorrect) is injected to the code for
presentation purposes.

To provide a more intuitive illustration, CEAM allows the programmer to query
AutoProof further to obtain a minimal counterexample in the last row of Fig.
2, where integer variables have been reduced to their minimal possible values.

Figure 2: Proof result in AutoProof: the first row (highlighted in red)
indicates a proof failure; the second row is a counterexample generated based
on the Z3 model; the third row is a minimized counterexample.

This minimized counterexample provides a simple diagnostic trace of max: on
loop initialization at line 7, Result = 0 and i = 2; at line 13, the exit
condition of the loop evaluates to True with a.count = 2 and i = 2, which
forces the loop to terminate. These values reveal the fault in the program:
the loop terminates too early, preventing the program from getting to the
actual maximum value, found at position 2 of the array a. To eliminate this
error, it suffices to strengthen the exit condition to permit one more loop
iteration: change i $\geq$ a.count to i $>$ a.count.

##  III Technology stack

This section introduces technologies used in the present work, including
language and prover.

Eiffel [15, 16] is an object-oriented programming language which natively
supports Design-by-Contract[14]. An Eiffel program consists of a set of
_classes_. A class represents a set of run-time objects characterized by the
features applicable to them. Fig. 3 shows a simple class representing bank
accounts. The class contains two types of features: _attributes_ representing
data items associated with instances of the class, such as balance (line 2)
and credit_limit (line 4); _routines_ representing operations applicable to
these instances, including available_amount and transfer. Routines are further
divided into _procedures_ (with no returned value) and _functions_ (returning
a value). Here, available_amount is a function returning an integer
(represented by the special variable Result), and transfer is a procedure.

⬇

1class ACCOUNT feature

2balance: INTEGER

3 \-- Balance of this account.

4credit_limit: INTEGER

5 \-- Credit limit of this account.

6available_amount: INTEGER

7 \-- Amount available on this account.

8 do

9 Result $:=$ balance - credit_limit

10 end

11transfer (amount: INTEGER; other: ACCOUNT)

12 \-- Transfer ‘amount’ to the ‘other’ account.

13 require

14 amount $>=$ $0$

15 amount $<=$ available_amount

16 do

17 balance $:=$ balance - amount

18 other.balance $:=$ other.balance $+$ amount

19 ensure

20 withdrawal: balance $=$ old balance - amount

21 deposit: other.balance $=$ old other.balance $+$ amount

22 end

23end

Figure 3: A class implementing the behavior of bank accounts

Programmers can specify the properties of Eiffel classes by equipping them
with contracts of the following types:

  * •

A precondition (require) must be satisfied at the time of any call to the
routine; the precondition of transfer (lines 13 – 15), for example, requires
the value of amount to be non-negative an no greater than available_amount.

  * •

A postcondition (ensure) must be guaranteed on routine’s exit; for instance, a
postcondition of transfer at line 20 states that, at the end of the execution
of transfer, the value of balance must have been decreased by amount.

  * •

A loop invariant (invariant) characterizes the semantics of a loop in the form
of a property satisfied after initialization and preserved by every iteration,
as illustrated by the invariant of max (lines 9 – 11 in Fig. 1) specifies the
properties of i and Result before and after every iteration.

  * •

A loop variant (variant) is an integer measure that should always be non-
negative and decrease strictly at each loop iteration, ensuring that the loop
eventually terminates; the loop variant of the loop in max is a.count - i
(line 19 in Fig. 1).

Contracts embedded in the code make it amenable to both dynamic analysis (run-
time checking of the contract properties), as in the EiffelStudio environment,
and static analysis (Hoare-style proofs), as in AutoProof.

AutoProof [1, 23] is a static verifier that checks the correctness of Eiffel
programs against their functional specifications (contracts).

When verifying an Eiffel program, AutoProof translates the program into a
Boogie program [2, 11], which is then transformed into a set of verification
conditions (VCs) in SMT-LIB [3] format, based on Dijkstra’s weakest
precondition calculus [8]. The program’s correctness follows from the validity
of the VCs. Boogie asks an SMT solver (by default Z3, as noted) to reason
about the validity of each VC. Specifically, the solver tries to find a model
(an interpretation of variables and functions used in the SMT encodings) that
satisfies the negation of a VC. If the solver is unable to find such a model
(no counterexample exists and thus VC is a tautology), the verification is
successful. If it succeeds in obtaining such a model, the verification fails
and the solver makes the model available. This model witnesses the invalidity
of the VC [12] and thus can be seen as a counterexample111The counterexample
is a potential counterexample since it can occasionally be spurious because of
the prover’s incompleteness, although that phenomenon is not significant in
our experience. at the SMT level. In general, an SMT model describes an
execution trace (a sequence of program states) of a failed routine, along
which the program goes to an error state. The CEAM approach makes use of such
models to generate easy-to-understand counterexamples.

##  IV Counterexample extraction and minimization

This section first shows how to generate counterexamples based on SMT models,
then presents the details of the CEAM strategy for counterexample
minimization.

###  IV-A Counterexample extraction

In general, to construct a counterexample it suffices to extract the concrete
values of relevant variables from the corresponding SMT model, and to use
these values to produce a counterexample message (as in Fig. 2). The format of
the message can vary depending on the chosen “verbosity level”. As the goal of
the approach and the tools is to ease the burden on programmers, the message
only displays the initial values of relevant input variables in the
counterexample, as illustrated in the case at the beginning of this article.

Fig. 4 shows a simplified portion of the Z3 model222As the models are too big
to include in their entirety, this presentation only displays the parts
relevant to the discussion. corresponding to the failed proof of transfer’s
postcondition labeled by _withdrawal_ (line 20 in Fig. 3). To construct a
counterexample for this failure, it suffices to obtain the initial values of
its three input variables: the implicit variable Current333Current represents
the active object in the current execution context, similar to this in Java.
and the two arguments amount and other.

⬇

1amount $\to$ $5799$

2Heap $\to$ T@U!val!$17$

3Current $\to$ T@U!val!$18$

4other $\to$ T@U!val!$18$

5ACCOUNT.balance $\to$ T@U!val!$7$

6ACCOUNT.credit_limit $\to$ T@U!val!$8$

7Select $\to$ {

8T@U!val!$17$ T@U!val!$18$ T@U!val!$7$ $\to$ $(-2147475928)$

9T@U!val!$17$ T@U!val!$18$ T@U!val!$8$ $\to$ $(-2147481727)$

10}

Figure 4: A slice of model of the proof failure of _withdrawal_

In the transformation from Eiffel program to SMT code, to encode the execution
semantics of an object-oriented program, the evolution of the _heap_ (the
collection of program objects) during an execution is modeled as a sequence of
constants prefixed with Heap. Here the SMT constant Heap (line 2) corresponds
to the heap at the initial program state of transfer. Select (line 7) is a
function for retrieving the values of objects’ fields. It takes three
parameters, i.e., a heap state, an object reference and a data field, and
returns the value of the specified field.

As the example shows, the concrete values of primitive variables (such as
amount) are given directly, whereas the values of non-primitive variables
(e.g., Current and other of ACCOUNT type) appear in a symbolic form, prefixed
with T@U!val!. Such symbolic values can be seen as abstract memory locations
for the corresponding variables (see [13, 4]). The Select function is
available to obtain the values of the corresponding fields. For example,
Current is an instance of ACCOUNT and thus has fields balance and
credit_limit. The initial value of balance can be obtained by applying Select
to a tuple made of Heap = T@U!val!17 (line 2), Current = T@U!val!18 (line 3)
and ACCOUNT.balance = T@U!val!7 (line 5); the tuple matches the mapping in
line 8, therefore the returned value for balance is $-2147475928$. Similarly,
the value of credit_limit can be retrieved through the mapping in line 9.

To display the value of a non-primitive variable in the resulting
counterexample message, the strategy first checks whether the variable has an
alias that has been looked up earlier; if so, it displays the alias relation
between the variable and its earliest alias in the message; otherwise, it
looks up all of its primitive data fields transitively and display them in the
message.

In this example, the model shows that other and Current have the same symbolic
values. Current is, consequently, an alias of other. As Current is processed
prior to other, the fields of other will not be looked up and the message will
display the alias relation between other and Current.

After applying the above rules, the counterexample can be derived: balance =
$-2147475928$, credit_limit = $-2147481727$, amount = $5799$, other = Current.

For variables of container types such as arrays or sequences, the resulting
message displays the values of their sizes and containing elements. Here, we
use the example of max in Fig. 1 to demonstrate counterexample extraction for
container types.

The AutoProof approach specifies container structures in terms of mathematical
structures [23]. For example, the content of the input array a of max is
specified through a special attribute sequence (see lines 10 – 11 in Fig. 1),
which represents the mathematical sequence of integer values stored in a’s
cells. To obtain the content of a from the counterexample, we need to get the
content of its sequence field from the model. Fig. 5 shows a slice of the
model for the proof failure of max. CEAM first extracts the value of sequence
by querying Select with the values of Heap, a, and ARRAY^INTEGER^.sequence
(lines 1 – 3). Those values match to the mapping in line 5, hence the value of
sequence is T@U!val!38. By using this value, CEAM then queries the two
functions Seq#Length and Seq#Item to get the values of sequence’s size (line
8) and elements (lines 12 – 14), respectively.

⬇

1Heap $\to$ T@U!val!$26$

2a $\to$ T@U!val!$18$

3ARRAY^INTEGER_32^.sequence $\to$ T@U!val!$9$

4Select $\to$ {

5T@U!val!$26$ T@U!val!$18$ T@U!val!$9$ $\to$ T@U!val!$38$

6}

7Seq#Length $\to$ {

8 T@U!val!$38$ $\to$ $11800$

9 T@U!val!$40$ $\to$ $28101$

10}

11Seq#Item $\to$ {

12 T@U!val!$38$ $1$ $\to$ $0$

13 T@U!val!$38$ $11799$ $\to$ $0$

14 T@U!val!$38$ $11800$ $\to$ $5$

15}

Figure 5: A snippet of the model of proof failure of _is_max_

###  IV-B Counterexample minimization

Some of the extracted values found to cause a failure, such as $-2147481727$
in the above counterexample, are too large to enable a programmer to visualize
easily the cause of the proof failure. CEAM can simplify counterexamples by
reducing the absolute value of such integers. Program verification is modular,
meaning it processes each routine independently; so does minimization.

As the counterexample of a routine $r$ consists of the initial states of its
input variables, to minimize a counterexample of $r$ it suffices to minimize
each of $r$’s input variables. The task of minimizing an input variable $x$
can be reduced to a set of integer minimization tasks based on the type of
$x$. The procedure minimize_integer in in Algorithm 1 minimizes an integer
variable. It applies to the absolute value, retaining the sign (lines 1 – 5).
If $x$ is an object reference, the algorithm first checks whether $x$ denotes
a container; if yes, it finds the minimal size of $x$ (lines 8 – 9) and then
minimizes $x$’s elements one by one (lines 11 – 13); otherwise it minimizes
each of $x$’s fields (lines 15 – 16).

1 if  _x is an integer_ then

2 if  _v $>$ 0_ then

3 minimize_integer ($x$)

4 elseif  _v $<$ 0_ then

5 minimize_integer ($-x$)

6

7 elseif  _x is an object reference_ then

8 if  _$x$ is a container_ then

9 minimize_integer ($x.count$)

10 $n\leftarrow$ minimized value of $x.count$

11 from

12 $i\leftarrow 1$

13 until

14 $i\leq n$

15 loop

16 minimize_general ($x[i]$)

17 $i\leftarrow i+1$

18

19 else

20 across each field $f$ of $x$ as $x.f$ loop minimize_general ($x.f$)

21

Algorithm 1 minimize_general ($x$): minimize a variable of integer or
reference type

Algorithm 2 shows the details of minimize_integer. $B$ represents the Boogie
procedure of routine $r$ generated by AutoProof. The core idea is to find the
smallest integer bound $m$ such that adding a precondition $0\leq x<m$ (line
10) to $B$ still yields the same verification results. When the algorithm ends
(no smaller value of $m$ can be found), the model from the last verification
run is the minimal possible.

1 from

2 $m$ $\leftarrow$ current value of $x$

3 $B$.add_precondtion ($0\leq x$ $\wedge$ $x<m$)

4 verify

5until

6 $no\ smaller\ value\ can\ be\ tried$

7loop

8 $B$.remove_last_precondition

9 $m$ $\leftarrow$ pick a smaller value

10 $B$.add_precondtion ($0\leq x$ $\wedge$ $x<m$)

11 verify

12

Algorithm 2 minimize_integer ($x$): minimize an integer

The algorithm starts by assigning to $m$ the value of $x$ in the initial
model, then iteratively decreases $m$. At each iteration, it adds a new
precondition to $B$ to reduce the the range of $x$ and performs verification
based on the updated $B$, to check whether there still exists a model with a
smaller value of $x$ within the interval [$0,\ m$).

Picking a smaller value for $m$ (line 9) can be implemented either by
sequentially decreasing the value or using a binary reduction (as in binary
search) for acceleration. The current implementation first checks whether $x$
can be 0; if yes, the minimization stops as it has found the minimum;
otherwise, it continues applying binary reduction. For more flexibility, it
uses two user-specified parameters controlling termination:

  * •

_Tolerance_ : lower bound on the size of interval used in the binary search
algorithm;

  * •

_Max iteration_ : maximum number of verification iterations allowed when
minimizing an integer.

We have not endeavored to prove the correctness of the algorithms since the
correctness of the approach (the “proof of the pudding”) is embodied in the
result: as the overall goal is to obtain a counterexample, the final criterion
is whether the minimized value is still a counterexample, as established
rigorously by the underlying proof technology. If not, the original
unminimized counterexample still applies.

##  V Experiment

A preliminary evaluation of the usability of the CEAM approach covers over 40
program versions444https://github.com/huangl223/Proof2Test/tree/main/examples
of 9 examples, including some adapted from software verification competitions
[24, 5, 10]. The examples (listed in Table I) include: 1) ACCOUNT introduced
in Fig. 3; 2) a CLOCK class implementing a clock counting seconds, minutes,
and hours; 3) a HEATER class implementing a heater adjusting it state (on or
off) based on the current temperature and a user-defined temperature; 4) a
LAMP class describing a lamp equipped with a switch (for switching on/off the
lamp) and a dimmer (for adjusting the light intensity of the lamp); 5)a
BINARY_SEARCH class implementing the binary search algorithm; 6) a
LINEAR_SEARCH implementing the linear search algorithm; 7) a SQUARE_ROOT that
calculates two approximate square roots of a positive integer; 8) MAX from
Fig. 1; 9) a SUM_AND_MAX class computing the maximum and sum of the elements
of an array.

Each row in the table reports on the experiment result of a single example,
which consists of multiple versions. Each version was intentionally injected
with a fault, such as confusions between $+$ and $-$, $\leq$ and $<$, $>$ and
$\geq$, missing loop invariant(s), pre- or postcondition, etc. The experiments
use AutoProof to verify the programs, produce counterexamples for all
occurring proof failures, and minimize them. The _tolerance_ and _max
iteration_ parameters are currently set to 0 and 20 respectively.

The third column gives the total number of integer variables whose minimized
in the experiment. Cases where no minimization is performed (e.g., the value
of an integer variable in the counterexample is already 0 before minimization)
are not included. The reduction rate (fourth column), number of iterations
(fifth column), verification time (sixth column) and minimization time (last
column) per integer are averaged out over all minimized integers of each
example.

As the experiment result shows, CEAM minimization is cost-effective overall:
in most cases, conducting minimization reduces the values of integer variables
by over 80% with an average cost of less than 4 extra verification runs
(iterations). Most of the minimized integer values are fairly small and easy-
to-understand: out of 125 minimized values, 108 are in the range [-2, 2], out
of which 58 are zero; values not in that range are usually close to the values
of some predefined constants in the program.

TABLE I: Experiment Results Example |  | Number
---
of versions
| Total Number of
---
Minimized Integers
| Avg. Reduction
---
Rate
| Avg. Number
---
of Iterations
| Avg. Verification
---
Time (seconds)
| Avg. Minimization
---
Time (seconds)
ACCOUNT | 7 | 17 | 99.98% | 2.5 | 0.028 | 0.087
CLOCK | 6 | 13 | 100% | 1.46 | 0.019 | 0.034
HEATER | 2 | 4 | 48.4% | 4.25 | 0.030 | 0.128
LAMP | 4 | 8 | 0.819% | 1.875 | 0.115 | 0.233
BINARY_SEARCH | 5 | 31 | 98.8% | 3.22 | 0.448 | 1.512
LINEAR_SEARCH | 3 | 9 | 99.9% | 3.44 | 0.087 | 0.279
SQUARE_ROOT | 4 | 3 | 89.9% | 4 | 0.133 | 0.505
MAX | 4 | 12 | 87.1% | 4.25 | 0.213 | 1.456
SUM_AND_MAX | 6 | 11 | 80.7% | 3.45 | 0.590 | 1.704

##  VI Related Work

In line with the objective of helping programmers to understand the causes of
proof failures, several approaches have been proposed to provide more user-
friendly visualizations of counterexample models [11, 6, 9, 22, 18]. Claire et
al. [11] developed the Boogie Verification Debugger (BVD), which interprets a
counterexample model as a static execution trace (i.e., a sequence of abstract
states). David et al. [9] transformed the models back into a counterexample
trace comprehensible at the original source code level (SPARK) and display the
trace using comments. Similarly, Aleksandar et al. [6] transformed SMT models
to a format close to the Dafny syntax. In contrast to the present work, these
approaches concentrate on the generation of human-readable counterexample and
do not consider any counterexample minimization.

Another direction of work to ease the understanding of proof failures is to
generate more useful counterexamples in the first place: Polikarpova et al.
[21] developed a tool, Boogaloo, which applies symbolic execution to generate
counterexamples for failed Boogie programs. Like the present approach,
Boogaloo displays minimal counterexamples in the form of valuations of
relevant variables. Müller et al. [17] implemented a Visual Studio dynamic
debugger plug-in for Spec#, to reproduce a failing execution from the view of
the prover. Likewise, Petiot et al. [20, 19] developed STADY, which produces
failing tests for the failed assertions using symbolic execution techniques.
That approach is also referred to as testing-based counterexample synthesis:
it first translates the original C program into programs suitable for testing
(run-time assertion checking), then applies symbolic execution to generate
counterexamples (input for failing tests) based on the translated program.
Unlike to that approach, CEAM counterexample extraction directly exploits the
counterexample models produced by the provers, and hence does not require
additional program instrumentation or counterexample generation.

##  VII Conclusion

This article has presented Counterexample Extraction and Minimization (CEAM),
an approach that improves the quality of counterexamples generated in the
presence of failed program proofs. CEAM automatically generates simple and
easy-to-understand counterexamples. We believe this makes the results of
failed proofs practical enough to be used by regular programmers. The CEAM
implementation is integrated in the AutoProof verifier to assist programmers
when debugging failed proofs. The approach could also be applied to other
Hoare-style verification tools relying on Boogie-style provers and SMT
solvers.

Ongoing work includes implementing a feature of automatic test generation
based on the counterexamples produced by CEAM, as well as extending the scope
of CEAM to include the supports for more data types and program constructs. We
also plan to conduct systematic empirical studies to evaluate precisely the
benefits of the proposed techniques for programmers with no verification
expertise.

Acknowledgments: We thank the anonymous referees for comments which led to
significant improvements. The work benefitted from discussions with Filipp
Mikoian, Alexander Kogtenkov and Alexander Naumchev from SIT.

## References

  * [1] AutoProof, http://comcom.csail.mit.edu/autoproof/
  * [2] Barnett, M., Chang, B.Y.E., DeLine, R., Jacobs, B., Leino, K.R.M.: Boogie: A Modular Reusable Verifier for Object-Oriented Programs. In: International Symposium on Formal Methods for Components and Objects. pp. 364–387. Springer (2005)
  * [3] Barrett, C., Stump, A., Tinelli, C., et al.: The SMT-LIB Standard: Version 2.0. In: International Workshop on Satisfiability Modulo Theories. vol. 13, p. 14 (2010)
  * [4] Bjørner, N., Moura, L.d., Nachmanson, L., Wintersteiger, C.M.: Programming Z3. In: International Summer School on Engineering Trustworthy Software Systems. pp. 148–201. Springer (2018)
  * [5] Bormer, T., Brockschmidt, M., Distefano, D., et al.: The COST IC0701 Verification Competition. In: International Conference on Formal Verification of Object-Oriented Software (FoVeOO). pp. 3–21. Springer (2011)
  * [6] Chakarov, A., Fedchin, A., Rakamarić, Z., Rungta, N.: Better Counterexamples for Dafny. In: International Conference on Tools and Algorithms for the Construction and Analysis of Systems (TACAS). pp. 404–411. Springer (2022)
  * [7] De Moura, L., Bjørner, N.: Z3: An Efficient SMT Solver. In: International conference on Tools and Algorithms for the Construction and Analysis of Systems (TACAS). pp. 337–340. Springer (2008)
  * [8] Dijkstra, E.W.: A Discipline of Programming. Prentice Hall (1976)
  * [9] Hauzar, D., Marché, C., Moy, Y.: Counterexamples from Proof Failures in SPARK. In: International Conference on Software Engineering and Formal Methods (SEFM). pp. 215–233. Springer (2016)
  * [10] Klebanov, V., Müller, P., , et al.: The 1st Verified Software Competition: Experience Report. In: International Symposium on Formal Methods (FM). pp. 154–168. Springer (2011)
  * [11] Le Goues, C., Leino, K.R.M., Moskal, M.: The Boogie Verification Debugger. In: International Conference on Software Engineering and Formal Methods. pp. 407–414. Springer (2011)
  * [12] Leino, K.R.M., Millstein, T., Saxe, J.B.: Generating Error Traces from Verification-Condition Counterexamples. Science of Computer Programming 55(1-3), 209–226 (2005)
  * [13] Leino, K.R.M., Rümmer, P.: The Boogie 2 Type System: Design and Verification Condition Generation
  * [14] Meyer, B.: Applying “Design by Contract”. Computer 25(10), 40–51 (1992)
  * [15] Meyer, B.: Object-Oriented Software Construction, vol. 2. Prentice Hall (1997)
  * [16] Meyer, B.: Touch of Class: Learning to Program Well with Objects and Contracts. Springer (2016)
  * [17] Müller, P., Ruskiewicz, J.N.: Using Debuggers to Understand Failed Verification Attempts. In: International Symposium on Formal Methods (FM). pp. 73–87. Springer (2011)
  * [18] Nilizadeh, A., Calvo, M., Leavens, G.T., Cok, D.R.: Generating Counterexamples in the Form of Unit Tests from Hoare-style Verification Attempts. In: International Conference on Formal Methods in Software Engineering (FormaliSE). pp. 124–128. IEEE (2022)
  * [19] Petiot, G., Kosmatov, N., Botella, B., Giorgetti, A., Julliand, J.: Your Proof Fails? Testing Helps to Find the Reason. In: International Conference on Tests and Proofs (TAP). pp. 130–150. Springer (2016)
  * [20] Petiot, G., Kosmatov, N., Botella, B., Giorgetti, A., Julliand, J.: How Testing Helps to Diagnose Proof Failures. Formal Aspects of Computing (FAC) 30(6), 629–657 (2018)
  * [21] Polikarpova, N., Furia, C.A., West, S.: To Run What No One Has Run Before: Executing an Intermediate Verification Language. In: International Conference on Runtime Verification (RV). pp. 251–268. Springer (2013)
  * [22] Stoll, C.: SMT Models for Verification Debugging. Master thesis, ETH Zurich (2019)
  * [23] Tschannen, J., Furia, C.A., Nordio, M., Polikarpova, N.: Autoproof: Auto-active Functional Verification of Object-Oriented Programs. In: International Conference on Tools and Algorithms for the Construction and Analysis of Systems (TACAS). pp. 566–580. Springer (2015)
  * [24] Weide, B.W., Sitaraman, M., Harton, H.K., Adcock, B., Bucci, P., Bronish, D., Heym, W.D., Kirschenbaum, J., Frazier, D.: Incremental Benchmarks for Software Verification Tools and Techniques. In: Working Conference on Verified Software: Theories, Tools, and Experiments (VSTTE). pp. 84–98. Springer (2008)


# A Second Order Thermal and Momentum Immersed Boundary Method for Conjugate
Heat Transfer in a Cartesian Finite Volume Solver

Ryan Crocker<EMAIL_ADDRESS>Yves Dubief  University of Vermont, Burlington,
Vermont, USA  Olivier Desjardins  Cornell University, Ithaca, New York, USA

###### Abstract

A conjugate heat transfer (CHT) immersed boundary (IB and CHTIB) method (IBM
and CHTIBM) is developed for use with laminar and turbulent flows with low to
moderate Reynolds numbers. The method is validated with the analytical
solution of the flow between two co-annular rotating cylinders at $Re=50$
which shows second order error convergence of the $L_{2}$ and $L_{\infty}$
error norms, of the temperature field, over a wide rage of solid to fluid
thermal conductivities, $\kappa_{s}/\kappa_{f}=\left(9-900\right)$. To
evaluate the CHTIBM with turbulent flow a fully developed, heated, turbulent
channel $\left(Re_{u_{\tau}}=150\text{ and }\kappa_{s}/\kappa_{f}=4\right)$ is
used which is correlated with to previous direct numerical simulation (DNS)
results. The CHTIBM is paired with a momentum IBM, both of which use a level
set field to define the wetted boundaries of the fluid/solid interfaces and
are applied to the flow solver implicitly with rescaling of the difference
operators of the finite volume (FV) method (FVM).

††journal: Journal of LaTeX Templates

##  1 Introduction

DNS has been a robust tool in numerical simulation and validation of new
algorithms for fluid flow and heat transfer by comparing experimental and DNS
results [1]. Much attention has been paid to canonical flows, such as the DNS
of fully developed turbulent channels in, [2] and [3], in these validation
studies. Fully developed flow and heat transfer have been studied
experimentally in [4], [5] and [6], and with DNS [7], [8], [9], and [10],
where numerical simulations were correlated with experimental data giving the
proper scaling and simulation parameters for turbulent, heated, flows.

The addition of CHT to DNS heat transfer studies greatly complicates the
analysis of developed flows and initial research into the subject suffered
from a lack of an adequate order of error. Such as in, [11] and [12], where
first order error convergence was observed. Recently, advances in the
interpolation schemes of curvilinear flow solvers has been shown to generate
solutions with second order error convergence, [13]. But this was shown for
small ratios of thermal conductivity and was not applied to the solution
implicitly or in a way easily applicable to cartesian FV solvers with complex,
non-mesh conforming boundary conditions.

Non-mesh conforming boundary conditions, IBMs, have been used to simulate
complex geometries with moderate Reynold’s number flows since their
development by Peskin,[14], to simulate heart valves. Yusof expanded the idea
of IBs and developed the direct forcing method, [15], eliminating velocity
stability issues associated with non-direct forcing IBMs. The final issues
with stability linked to small cut cells and pressure oscillations [16] were
solved with by the reconstruction or cut-cell IBMs where multiple methods were
developed: cell merging [17], mixed cell linking/merging [18], cell
linking[19], and cell mixing [20] methods, all of which show second order
error in their results for low to moderate Reynold’s numbers. In this study
the sharp representation of boundaries is done through a mixing method
modified from [20].

In the present work a DNS flow solver, developed in [21], with CHTIB and
momentum IB algorithms, with scalar transport, is used to study the effect of
CHT on two separate canonical flows. With the domains that are separated with
level set fields from which the cut-cell geometric data is derived. The first
being a heated co-annular cylinder set up in two dimensions with a simple
analytical solution which shows second order accuracy for very large
solid/fluid conductivity ratios at low Reynolds numbers. The second half of
this study uses a fully developed channel heated from the top and the bottom
boundaries where the root-mean-square of the near-boundary fluctuating
temperature is shown to correlate well with previous results from [22].

##  2 Governing Equations

The basic equations governing the transport of momentum and the scalar
(temperature) carried by the flow are given in their most general forms in
sections 2.1 and 2.2. The are kept fully conservative and will subsequently be
non-dimensionalized for each numerical experiment.

###  2.1 Momentum and mass conservation

The two equations used to simulate fluid flow are the Navier-Stokes set of
equations,

$\partial_{t}\left(\rho\textbf{u}\right)+\boldsymbol{\nabla}\cdot\rho\left(\textbf{u}\otimes\textbf{u}\right)=-\boldsymbol{\nabla}\cdot\textbf{I}P+\boldsymbol{\nabla}\cdot\left(\mu\boldsymbol{\nabla}\textbf{u}\right),$
(1)

and the conservation of mass, or continuity equation,

$\partial_{t}\left(\rho\textbf{u}\right)+\boldsymbol{\nabla}\cdot\left(\rho\textbf{u}\right)=0.$
(2)

In equations 2 and 1, bold terms are vectors, u is the velocity vector
$\left(\textbf{u}=u\hat{\boldsymbol{\imath}}+v\hat{\boldsymbol{\jmath}}+w\hat{\textbf{k}}\right)$,
$P$ is the pressure, $\mu$ is the dynamic viscosity, $\rho$ is the density, I
is identity matrix, and $\partial_{t}$ is the temporal derivative.

###  2.2 Scalar transport

Temperature is transported as a passive scaler in the flow, via equation 3,

$\partial_{t}\phi+\boldsymbol{\nabla}\cdot\left(\textbf{u}\phi\right)=\boldsymbol{\nabla}\cdot\left(\alpha_{f}\boldsymbol{\nabla}\phi\right),$
(3)

where $\phi$ is the temperature, and $\alpha$ is the thermal diffusivity. For
heat transfer inside the solid domain the convective term is dropped from
equation 3 giving,

$\partial_{t}\phi=\boldsymbol{\nabla}\cdot\left(\alpha_{s}\boldsymbol{\nabla}\phi\right),$
(4)

where the $f$ and $s$ subscripts stand for the fluid and solid domains,
respectively.

##  3 Sub-Domain by Level Set Field

The fluid and solid domains are separated by different signed regions of a
level set field [23], $G$, over the zero isosurface, $\Gamma\rightarrow G=0$.
The use of a particular form of the level set field, a signed distance
function, allows the immersed boundaries to represent a non-mesh conforming
boundary on a cartesian grid. The main property of a signed distance function,
$\boldsymbol{\nabla}G|=1$, allows this by giving a definite a zero iso-
surface, $G=0\rightarrow\Gamma$, where at all times a normal is defined by
$\textbf{n}=\nabla G/|\nabla G|$, and easily accessible geometric information
for cut-cell areas (see figure 2(b)) and volumes needed in section 4. The
geometric information is found with a simplex marching algorithm which uses
the intersection of tetrahedra with $\Gamma$. This information is then used to
locate the cut cell intersections with $\Gamma$ and the resulting cut-cell
geometry [24]. The last caveat of level set domain definition is that sign of
the $G$ also sets the parameters of each subdomain, where $G\geq 0$ is the
solid subdomain and $G<0$ is the fluid domain, but the choice of this is
arbitrary.

##  4 Numerical Methods

The general conservation form of a transport equation is a typical first step
in the derivation of the finite volume method. As the IB methods for both
energy and momentum are themselves derived from the finite volume method
equations 2 and 3 are rearranged into general conservation equations on a
staggered pressure/velocity grid (see figures 1(a) and 2(a)). This begins with
definition of the flux function, $F$, for momentum the flux function is a
tensor,

$[\textbf{F}\left(P,\textbf{u}\right)]=\left(\rho\textbf{u}\otimes\textbf{u}\right)+\textbf{I}P-\left(\mu\boldsymbol{\nabla}\textbf{u}\right),$
(5)

and for temperature transport it is a vector,

$\textbf{F}\left(\phi,\textbf{u}\right)=\left(\alpha\boldsymbol{\nabla}\phi\right)-\left(\textbf{u}\phi\right),$
(6)

Equations 5 and 6 are now in the form of general conservation equations,

$\underbrace{\partial_{t}\Theta+\boldsymbol{\nabla}\cdot\textbf{F}=0}_{energy}\hskip
21.33955pt\text{or}\hskip
21.33955pt\underbrace{\partial_{t}\boldsymbol{\Theta}+\boldsymbol{\nabla}\cdot[\textbf{F}]=0}_{momentum}$
(7)

The equations in 7 are also known as the advective form of the linear
conservation law, balance between the diffusive transport budget and the
convective transport budget, for time $t$, gives the overall change of the
conserved quantity, $\Theta$, as a function of time.

###  4.1 Basic finite volume discretization

Using the fully conservative form of the momentum equations gives a second
order accurate discretization of equation 1 , as shown in [25], [26] and [21],
for low mach number flows. Discretizing equation 3 in the same manner as that
shown in [27] will also give second order accuracy, with high-ordered-
upwinded-centered scheme (HOUC), [28]. The basic finite volume (FV) method
(FVM) can most generally be derived with the momentum form of the conservation
equation given in equation 7.

(a)

(b)

Figure 1: a.) A schematic of a 2D cut cell with the IB surface shown. b.) A
schematic of the small cell mixing used to fix the cells that do not fall
under the $\epsilon$ criteria.

Integrating equation 7 over the C.V. in figure 2(a) and replacing the
conserved quantity, $\boldsymbol{\Theta}$ with u and dropping the tensor
brackets around $F$, gives the integral form of the conservation law for
momentum,

$\int_{V}\partial_{t}\left(\rho\textbf{u}\right)dV+\int_{V}\boldsymbol{\nabla}\cdot\textbf{F}dV=0\hskip
7.11317pt\rightarrow\hskip
7.11317pt\partial_{t}\overline{\textbf{u}}=\frac{1}{V}\oint_{S}\textbf{F}\cdot\textbf{n}dA,$
(8)

where $V=\Delta x\Delta y\Delta z$,
$\overline{\textbf{u}}=\int_{V}\textbf{u}dV$ is the average velocity over the
volume, $A$ is the side area and n the cell face normal from the C.V. in
figure 2(a), and Gauss’s law has been used from left to right. Omitting the
over bar on u, and integrating around the surface $S$, the expanded, spatially
discretized form of equation 8 is,

$\begin{split}\partial_{t}\left(\rho\textbf{u}\right)=\frac{1}{V}\bigg{[}&\left(\textbf{Fx}_{i+\frac{1}{2},j,k}A_{i+\frac{1}{2},j,k}-\textbf{Fx}_{i-\frac{1}{2},j,k}A_{i-\frac{1}{2},j,k}\right)\cdot\textbf{n}_{x}\\\
+&\left(\textbf{Fy}_{i,j+\frac{1}{2},k}A_{i,j+\frac{1}{2},k}-\textbf{Fy}_{i,j-\frac{1}{2},k}A_{i,j-\frac{1}{2},k}\right)\cdot\textbf{n}_{y}\\\
+&\left(\textbf{Fz}_{i,j,k+\frac{1}{2}}A_{i,j,k+\frac{1}{2}}-\textbf{Fz}_{i,j,k-\frac{1}{2}}A_{i,j,k-\frac{1}{2}}\right)\cdot\textbf{n}_{z}\bigg{]},\end{split}$
(9)

where $A$ is the area of the cell face, $A_{i,j\pm\frac{1}{2},k}=\Delta
x\Delta z$, $A_{i\pm\frac{1}{2},j,k}=\Delta y\Delta z$,
$A_{i,j,k\pm\frac{1}{2}}=\Delta x\Delta y$, $\textbf{n}=[\textbf{n}_{x}\hskip
7.11317pt\textbf{n}_{y}\hskip 7.11317pt\textbf{n}_{z}]$ are the normal
components from each face, and $\textbf{F}=\left[\textbf{Fx}\hskip
7.11317pt\textbf{Fy}\hskip 7.11317pt\textbf{Fy}\right]$ are the flux functions
components in the $x$, $y$, and $z$ $\left(\hat{\boldsymbol{\imath}}\hskip
7.11317pt\hat{\boldsymbol{\jmath}}\hskip 7.11317pt\hat{\textbf{k}}\right)$
directions.

(a)

(b)

Figure 2: a.) A computational, finite volume, cell for momentum, with
temperature and pressure carried at the center node. b.) A 3D cut cell for the
thermal or momentum IB showing a more general 3D cut cell.

Replacing the flux function and the time derivative in equation 9 with the
flux function from equation 6 and the time derivative from equation 3,
respectively, gives the FV formulation of the energy transport equation.

###  4.2 Conservative Temporal Discretization

Time advancement of equation 9 is done by use of a 3 directional semi-implicit
fractional step method, [29], using a second order in time Crank-Nicholson
discretization scheme, with approximate factorization, [30], which is a form
of alternating direction implicit methods (ADIM). The ADIM system of equations
is trilinear and solution of which is readily expanded to three dimensions,
and to codes written in parallel with message passing interfaces (MPI),

$\partial\left(\rho\textbf{u}\right)=\frac{\left(\rho\textbf{u}\right)^{n+1}-\left(\rho\textbf{u}\right)^{n}}{\Delta
t}=\textbf{L}\left(\textbf{u},P\right),$ (10)

where $n$ denotes the present time step, $n+1=t+\Delta t$, and L is the linear
operator
$\textbf{L}\left(\textbf{u},P\right)=-\boldsymbol{\nabla}\cdot\textbf{F}\left[\frac{1}{2}\left(\textbf{u}^{n+1}+\textbf{u}^{n}\right),P^{n+1}\right]$.
Newton-Iterations, $m$, are applied to the sub-steps of equation 10, for each
time step, $\Delta t$. The method is taken from [31] and is used here, as an
example in the $x$ direction, as,

$\begin{split}&\left[\textbf{I}-\frac{1}{2}\Delta
t\partial_{\textbf{u}}\textbf{L}\left(\textbf{u},P\right)\right]\left(\rho
u^{n+1}_{m+1}-\rho u^{n}_{m}\right)=\\\ &\rho u^{n}_{m=0}-\rho
u^{n}_{m}+\Delta
t\textbf{L}\left[\frac{1}{2}\left(u^{n}+u^{n+1}\right),Px^{n+1}_{m}\right],\end{split}$
(11)

where $m=0$ is the initial solution at the present time step, $n$.

The velocity field resulting from the solution of equations 10-11, using
equation 10 9, will not be solenoidal. To get a divergence free velocity
equation 2 is used to project the pressure onto a divergence free velocity
field give the Pressure-Poisson equation,

$\nabla^{2}\eta=-\boldsymbol{\nabla}\cdot\textbf{u}^{*}\rightarrow\textbf{u}^{n+1}=\textbf{u}^{n}-\nabla\eta\hskip
7.11317pt\text{and}\hskip 7.11317ptP^{n+1}=P^{n}+\frac{\eta}{\Delta t}$ (12)

where $u*=u^{n+1}$ is the intermediate velocity from the fractional step
method, $P^{n+1}$ is the updated pressure, and $\eta$ is the pseudo-pressure
field also from the fractional step method.

The time advancement scheme given in equations 10 and 11 is directly
applicable to the scalar transport spatial discretization given given in
sections 4 and 4.1. With the only difference being the replacement of the flux
function in the linear operator L with
$\textbf{L}\left(\phi,\textbf{u}\right)=-\boldsymbol{\nabla}\cdot\textbf{F}\left(\phi,\textbf{u}\right)$,
with the flux function from equation 6, and the fractional step method is not
used, as the continuity equation is not applied.

###  4.3 Momentum Immersed Boundary Method

The second order accurate, momentum cut-cell IB will be covered, briefly, in
this section. Much more in depth coverage of the topic can be found in [20],
[24] and [32]. Nominally, the IB methods for momentum and thermal transport
provide the same functionality. They add cells cut by a non-grid-conforming
boundary condition back into the time evolution solutions, minimizing error.
Where they differ in that the thermal IB from section 4.4 requires the
information on both sides of the surface $\Gamma$, where the momentum cut-cell
IB, given here, is only concerned with the wetted side of $\Gamma$.

Cutting the control volume cells adds two new terms to the momentum equations,
the viscous, $D$, and momentum, $C$, source terms. $D$ and $C$ enforce the no-
slip condition tangental to $\Gamma$, and the Lagrangian velocity of $\Gamma$,
respectively. To fully describe the FVM with cut cells, along with $D$ and
$C$, the cut cell aperture size, $\left(\zeta A\right)_{i,j,k}$, the area
fraction $\zeta_{i,j,k}$, and the cut cell volume fraction, $\lambda_{i,j,k}$,
and $V_{i,j,k}=\lambda_{i,j,k}\Delta x\Delta y\Delta z$ is the cell volume.
The addition of the IB will add the two new terms to equation 9,

$\begin{split}\partial_{t}\left(\rho\textbf{u}\right)=&+\frac{1}{\lambda_{i,j,k}\Delta
x}\left[\left(\zeta
A\right)_{i+\frac{1}{2},j,k}\textbf{Fx}^{n}_{i+\frac{1}{2},j,k}-\left(\zeta
A\right)_{i-\frac{1}{2},j,k}\textbf{Fx}^{n}_{i-\frac{1}{2},j,k}\right]\cdot\textbf{n}_{x}\\\
&+\frac{1}{\lambda_{i,j,k}\Delta y}\left[\left(\zeta
A\right)_{i,j+\frac{1}{2},k}\textbf{Fy}^{n}_{i,j+\frac{1}{2},k}-\left(\zeta
A\right)_{i,j-\frac{1}{2},k}\textbf{Fy}^{n}_{i,j-\frac{1}{2},k}\right]\cdot\textbf{n}_{y}\\\
&+\frac{1}{\lambda_{i,j,k}\Delta z}\left[\left(\zeta
A\right)_{i,j,k+\frac{1}{2}}\textbf{Fz}^{n}_{i,j,k+\frac{1}{2}}-\left(\zeta
A\right)_{i,j,k-\frac{1}{2}}\textbf{Fz}^{n}_{i,j,k-\frac{1}{2}}\right]\cdot\textbf{n}_{y}\\\
&+\frac{1}{\lambda_{i,j,k}V_{i,j,k}}\left[C+D\right].\end{split}$ (13)

Both $\zeta_{i,j,k}$ and $\lambda_{i,j,k}$ are $0<\zeta_{i,j,k}\text{ ,
}\lambda_{i,j,k}\leq 1$ and are the fraction of the cell face and volume lost
due to the cell cut (see figures 1(a) and 2(b)).

The viscous term added by the IB method, $D$, comes from the integration of
the equation 8 taking into account the interface $\Gamma$ from figure 1(a),
giving,

$\ \int_{\Gamma}=\nu
D=\nu\left(\nabla\textbf{u}\right)\cdot\textbf{n}_{\Gamma}dA_{\Gamma}.$ (14)

Equation 14 accounts for the no slip boundary condition at the surface cut by
the boundary. The momentum exchange term, $C$, accounts for any motion of the
interface cut. If the boundary is moving $C$ has the form,

$\ C=v_{\textbf{n}}\frac{V_{i,j,k}}{\Delta t}.$ (15)

If the interface boundary is not moving and impervious, $C=0$. When the
boundary is in motion the wall normal velocity, $v_{\textbf{n}}$ is non-zero.

The application of the cut cell method can lead to very small cell volumes
compared to the rest of the mesh. By not taking account of these small cells
the velocity and pressure solutions will be numerically unstable without,
small, and computationally prohibitive, time steps. The most efficient and
stable numerical treatment of the small cells is by merging, [33], them with
the surrounding cut cells and un-cut cells. Though giving a numerically stable
solution cell merging does not lend itself to straight forward implementation
with IB’s that are not stationary. A more versatile method is to mix the small
cell values with the cells that surround them such as the method used in [34].
The mixing method is altered to fit a staggered grid and again follows the
method in [20].

The mixing fraction, $\beta^{i}_{i,j,k}$ and the exchange variable, $\xi_{i}$,
are given in the $x$ direction with the conservative mixing variable,
$\Upsilon$, in equation 16

$\
\xi_{x}=\frac{\beta^{x}_{i,j,k}}{\beta^{x}_{i,j,k}V_{i,j,k}+V_{\text{target}_{x}}}\left(V_{i,j,k}\Upsilon^{*}_{\text{target}_{x}}-V_{\text{target}_{x}}\Upsilon^{*}_{\text{target}_{x}}\right).$
(16)

Mixing fractions, $\beta^{M}_{i,j,k}$, are a combination of the normal vector
to the target cell and $\lambda_{i,j,k}$, which in two dimensions is,

$\begin{split}&\beta^{M}_{i,j,k}=\textbf{n}^{2}_{M}\lambda_{i,j,k}\rightarrow\\\
&\beta^{x}_{i,j,k}=\textbf{n}^{2}_{\Gamma
x}\lambda^{\text{target}_{\text{\scalebox{1.5}{$x$}}}}_{i,j,k}\hskip
14.22636pt\beta^{y}_{i,j,k}=\textbf{n}^{2}_{\Gamma
y}\lambda^{\text{target}_{\text{\scalebox{1.5}{$y$}}}}_{i,j,k}\hskip
14.22636pt\beta^{xy}_{i,j,k}=\textbf{n}_{\Gamma x}\textbf{n}_{\Gamma
y}\lambda^{\text{target}_{\text{\scalebox{1.5}{$xy$}}}}_{i,j,k}.\end{split}$
(17)

The cut off for cells to be considered small for this work is
$\lambda_{i,j,k}\leq 0.1$, and the condition that the sum of all of the mixing
fractions be normalized so that, $\sum^{N}_{M=1}\beta^{M}_{i,j,k}=1$, where
$N$ is the number of target cells. The values for $\Upsilon$ are mixed to the
surrounding cells from the small cells by equation 18,

$\ \Upsilon=\Upsilon^{*}+\frac{1}{V_{i,j,k}}\sum^{N}_{M=1}\xi_{M},$ (18)

where in this case $\Upsilon^{*}$ is the variable to be mixed to before
mixing, and $N$ is the number of surrounding cells to the small cell. The form
of $\Upsilon^{*}$ that is gives the best solution to the pressure is,

$\ \Upsilon^{*}=\frac{1}{V_{i,j,k}}\int_{V_{i,j,k}}\textbf{F}\cdot\textbf{n},$
(19)

and will also be the method in which the small cell conservative fluxes are
mixed to their neighbors. After mixing, the fully discretized Pressure-Poisson
equation (equation 12), from the fractional step method, is written as
equation 20,

$\begin{split}&\frac{1}{\Delta x}\left[\left(\zeta
A\right)_{i+\frac{1}{2},j}\frac{\partial\eta}{\partial
x}\Bigg{|}_{i+\frac{1}{2},j}-\left(\zeta
A\right)_{i-\frac{1}{2},j}\frac{\partial\eta}{\partial
x}\Bigg{|}_{i-\frac{1}{2},j}\right]\\\ +&\frac{1}{\Delta y}\left[\left(\zeta
A\right)_{i,j+\frac{1}{2}}\frac{\partial\eta}{\partial
y}\Bigg{|}_{i,j+\frac{1}{2}}-\left(\zeta
A\right)_{i,j-\frac{1}{2}}\frac{\partial\eta}{\partial
y}\Bigg{|}_{i,j-\frac{1}{2}}\right]-\frac{\lambda_{i,j}\Gamma_{i,j}\nabla\eta\cdot\textbf{n}}{V_{ij}}\\\
=-&\frac{1}{\Delta x}\left[\left(\zeta
A\right)_{i+\frac{1}{2},j}u^{*}_{i+\frac{1}{2},j}-\left(\zeta
A\right)_{i-\frac{1}{2},j}u^{*}_{i-\frac{1}{2},j}\right]\\\ -&\frac{1}{\Delta
y}\left[\left(\zeta
A\right)_{i,j+\frac{1}{2}}v^{*}_{i,j+\frac{1}{2}}-\left(\zeta
A\right)_{i,j-\frac{1}{2}}v^{*}_{i,j-\frac{1}{2}}\right]+\frac{\lambda_{ij}\Gamma_{i,j}\text{v}_{\textbf{n}_{x}}}{V_{i,j}},\end{split}$
(20)

where this equation is easily expanded to three dimensions.

###  4.4 Thermal Immersed Boundary Method

The CHTIB derived form the conjugate heat transfer boundary conditions at the
interface of two dissimilar substances without contact resistance. Those
boundary conditions are the continuity of the normal surface flux rates at
$\Gamma$, shown in variable and discretized form,

$\kappa_{s}\frac{\partial\phi}{\partial\textbf{n}_{\Gamma
s}}=\kappa_{f}\frac{\partial\phi}{\partial\textbf{n}_{\Gamma
f}}\rightarrow\kappa_{s}\frac{3\phi_{\Gamma
s}-4\tilde{\phi}_{s1}+\tilde{\phi}_{s2}}{\partial\textbf{n}_{\Gamma
s}}=\kappa_{f}\frac{-3\phi_{\Gamma
f}+4\tilde{\phi}_{f1}-\tilde{\phi}_{f2}}{\partial\textbf{n}_{\Gamma f}},$ (21)

respectively, and the continuity of temperature at $\Gamma$,

$\phi_{\Gamma s}=\phi_{\Gamma f}=\phi_{\Gamma},$ (22)

where the $\tilde{\left(\cdot\right)}$ terms are the interpolated equidistant
(spaced at $\gamma_{\text{int}}$) from the interface on either side of the
domain using second order interpolation scheme (see figure 3(a)) s, and
$\kappa$ is the thermal conductivity. The discretized forms of the CHT
boundary conditions, given in equations 21 and 22, are used to form an
implicit, single boundary condition for any cut cell. Substituting equation
22, and the relationship $\textbf{n}_{\Gamma s}=-\textbf{n}_{\Gamma
f}=\textbf{n}_{\Gamma}$, into equation 21 gives the new flux condition,

$\left.\kappa_{s}\frac{3\phi_{\Gamma}-4\tilde{\phi}_{s1}+\tilde{\phi}_{s2}}{\partial\textbf{n}_{\Gamma}}\right|_{s}=\left.\kappa_{f}\frac{-3\phi_{\Gamma}+4\tilde{\phi}_{f1}-\tilde{\phi}_{f2}}{\partial\textbf{n}_{\Gamma}}\right|_{f}.$
(23)

Equation 23 is solved for $\phi_{\Gamma}$ which is used to apply the CHT
boundary condition to the cell containing $\phi_{i}$ from figure 3(b).

(a)

(b)

Figure 3: a.) One dimensional interpolation stencils used to find
$\phi_{\Gamma}$. b.) One dimensional stencil used to apply the boundary
condition from $\Gamma$ to $\phi_{i}$.

The CHTIB is implemented in the same manner as embedded FV boundary
conditions. Solving for $\phi_{\Gamma}$, using the interface distance,
$\gamma$ from figure 3(b), the Crank-Nicholson, semi-implicit form, of the
scalar transport equation,

$\displaystyle\partial_{t}\phi=\boldsymbol{\nabla}\cdot\left(\boldsymbol{\nabla}\alpha\phi\right)-\boldsymbol{\nabla}\cdot\left(\textbf{u}\phi\right)\rightarrow$
$\displaystyle
a_{i}\phi^{n+1}_{i-1}+b_{i}\phi^{n+1}_{i}+c_{i}\phi^{n+1}_{i+1}=a_{i}\phi^{n}_{i-1}+b_{i}\phi^{n}_{i}+c_{i}\phi^{n}_{i+1},$
(24)

where $a$, $b$, and $c$ are the linear operator metrics from the FV
discretization showing both the implicit, left hand side (LHS), and explicit,
right hand side (RHS). The embedded, one dimensional, boundary condition for
figure 3(b),

$\phi^{n}_{i+1}=\phi^{n}_{\Gamma}\frac{\Delta
x_{i}}{\gamma}+\phi^{n}_{i}\left(1-\frac{\Delta x_{i}}{\gamma}\right)$ (25)

given the implicit boundary conditions as,

$a_{i}\phi^{n+1}_{i-1}+\phi^{n+1}_{i}\left[b_{i}-\left(1-\frac{\Delta
x_{i}}{\gamma}\right)\right]\rightarrow c_{i}=0,$ (26)

and explicit as equation 27,

$\boldsymbol{\nabla}\cdot\left(\boldsymbol{\nabla}\alpha\phi\right)-\boldsymbol{\nabla}\cdot\left(\textbf{u}\phi\right)\hskip
14.22636pt\text{with}\hskip
14.22636pt\phi^{n}_{i+1}=\phi^{n}_{\Gamma}\frac{\Delta
x_{i}}{\gamma}+\phi^{n}_{i}\left(1-\frac{\Delta x_{i}}{\gamma}\right),$ (27)

in the flux terms. By replacing $\gamma$ in the boundary condition
interpolation, equation 25 it can be noted that equations 26 and 27 reduced to
standard, second order embedded boundary conditions using ghost nodes. It
should also be noted that the values of either right or left boundary
conditions, $\phi_{i+1}$ or $\phi_{i-1}$ are not replaced in the solution,
only their flux terms are affected.

There are some stability issues with the final form of equations 26 and 27.
The first, is $\lim_{\gamma\rightarrow 0}\frac{\Delta}{\gamma}=\infty$ causes
solution instability when $\gamma/\Delta x\leq 0.05$. By keeping the equations
26 and 27 in their given form for $\gamma/\Delta x\geq 0.05$, and replacing
the interpolation of the embedded boundary condition in equation 25 with
$\phi_{i}=\phi_{\Gamma}$, removes this stability constraint by adding a simple
Dirichlet boundary condition to the point collocated on $\Gamma$.

The second stability concern is how to treat $\phi^{n}_{\Gamma}$ and
$\phi^{n+1}_{\Gamma}$. It is clear from the form of equation 26 that it has
already been taken that $\phi^{n}_{\Gamma}=\phi^{n+1}_{\Gamma}$. This was
found to be consistent with analytical solutions, and using a $dt$ of the
order used for the smallest scales of $\phi$, will give second order results
(see section 5).

The third, and final stability constraint has to do with implementation in $N$
dimensions. When a FV cell is cut only the fluxes in the direction that will
cross $\Gamma$ is used to determine the form of the evolution equation for
that cell. For example, in the cells in shown in figure 4, in figure 4(a) the
CHT boundary condition would be applied for the $\phi_{i+1,j}$ and the fluxes
in the $y$ would be set equal at each face of the FV cell at $\\{i,j-1\\}$ and
$\\{i,j+1\\}$. In figure 4(b) the boundary lies across both that FV cell faces
connect $\phi_{i,j}$ to $\phi_{i+1,j}$ and $\phi_{i,j-1}$. The CHT boundary
conditions would then be applied at both of those points. Canceling out the
fluxes, near $\Gamma$, that do not apply directly to a CHT boundary for that
cell is justified by the physical assumptions that: near the interface the
fluxes not orthogonal to $\ \textbf{n}_{\Gamma}$ will not conduct heat across
that interface,

$\textbf{F}\times\textbf{n}_{\Gamma}=0\hskip 7.11317pt\text{if}\hskip
7.11317pt\textbf{F}\parallel\textbf{n}_{\Gamma},$ (28)

constraining the near interface physics to CHT. Second, iso-therms near
$\Gamma$ will be parallel to it giving,

$\boldsymbol{\nabla}\phi\parallel\textbf{n}_{\Gamma}<<\boldsymbol{\nabla}\phi\perp\textbf{n}_{\Gamma}.$
(29)

Thirdly, the simple numerical constraint that there is no sub-grid scale model
being used for clipped cells, other than the CHT one laid previously.

(a)

(b)

Figure 4: a.) Two dimensional schematic showing how $\Gamma$ would forced the
CHT boundary conditions to be applied at $\phi_{i+t,j}$. b.) Two dimensional
schematic showing how $\Gamma$ would forced the CHT boundary conditions to be
applied at $\phi_{i+t,j}$ and $\phi_{i,j-1}$.

Essentially the CHTIB decouples the fluid and solid domains by leveraging the
the FVM. The resulting form transfers information with high order
interpolation methods which pass through the CHT boundary conditions. In this
way, either domain can only indirectly affect the other, and only in a way
that is prescribed in equations 21-22.

##  5 Thermal IB Results and Validation

###  5.1 Co-Annular Validation

Validation was done with one of the few CHT analytical solution, using the
normalization of equations 1-2, for incompressible flow, given as,

$\partial_{t^{*}}\left(\textbf{u}^{*}\right)+\boldsymbol{\nabla}^{*}\cdot\left(\textbf{u}^{*}\textbf{u}^{*}\right)=-\boldsymbol{\nabla}^{*}\cdot\textbf{I}P^{*}+\frac{1}{Re}\boldsymbol{\nabla}^{*}\cdot\left(\boldsymbol{\nabla}^{*}\textbf{u}^{*}\right),$
(30)
$\partial_{t^{*}}\left(\textbf{u}^{*}\right)+\boldsymbol{\nabla}^{*}\cdot\left(\textbf{u}^{*}\right)=0.$
(31)

The equations for heat transfer, equations 3-4, are normalized as,

$\partial_{t^{*}}\phi^{*}+\boldsymbol{\nabla}^{*}\cdot\left(\textbf{u}^{*}\phi^{*}\right)=\frac{1}{RePr}\boldsymbol{\nabla}^{*}\cdot\left(\boldsymbol{\nabla}^{*}\phi^{*}\right),$
(32)
$\partial_{t}\phi^{*}=\frac{1}{RePr}\frac{\alpha_{s}}{\alpha_{f}}\boldsymbol{\nabla}^{*}\cdot\left(\boldsymbol{\nabla}^{*}\phi^{*}\right),$
(33)

where $Re$ is the Reynolds number and $Pr$ is the Prandtl number. The
dimensionless groups used in the normalizations in equations 30-31 and 32-33
are,

$\begin{split}\textbf{u}^{*}=\frac{\textbf{u}}{U_{\infty}}\end{split}\hskip
14.22636pt\begin{split}\textbf{x}^{*}=\frac{\textbf{x}}{D_{o}}\end{split}\hskip
14.22636pt\begin{split}t^{*}=t\frac{U_{\infty}}{D_{o}}\end{split}\hskip
14.22636pt\begin{split}P^{*}=\frac{P}{\rho U_{\infty}^{2}}\end{split}\hskip
14.22636pt\begin{split}Re=\frac{U_{\infty}D_{o}}{\nu},\end{split}$ (34)

where $D_{o}=2R_{o}$, and $U_{\infty}=u_{r}\left(R_{o}\right)$, defined in
figure 5 and equations 35 and 36.

The velocity and temperature fields for two co-annular cylinders (see figure
5(a)) are given by,

$u_{\theta}=0,\hskip 7.11317pt\text{and}\hskip
7.11317ptu_{r}\left(r\right)=\begin{dcases}0&\text{for}\hskip
5.69046ptR_{i}<r<R_{m}\\\
-\frac{R_{o}R^{2}_{m}U_{\infty}}{R^{2}_{o}-R^{2}_{m}}\frac{1}{r}+\frac{R_{o}U_{\infty}}{R^{2}_{o}-R^{2}_{m}}r&\text{for}\hskip
5.69046ptR_{m}<r<R_{o}\end{dcases}$ (35)

in cylindrical coordinates, and,

$\phi_{r}\left(r\right)=\begin{dcases}\phi_{in}+\frac{\phi_{out}-\phi_{in}}{\text{log}_{10}\left(\frac{R_{m}}{R_{i}}\right)+\left(\frac{\kappa_{s}}{\kappa_{f}}\right)\text{log}_{10}\left(\frac{R_{o}}{R_{m}}\right)}\text{log}_{10}\left(\frac{r}{R_{i}}\right)&\text{for}\hskip
5.69046ptR_{i}<r<R_{m}\\\
\phi_{out}-\frac{\phi_{out}-\phi_{in}}{\left(\frac{\kappa_{f}}{\kappa_{s}}\right)\text{log}_{10}\left(\frac{R_{m}}{R_{i}}\right)+\text{log}_{10}\left(\frac{R_{o}}{R_{m}}\right)}\text{log}_{10}\left(\frac{R_{o}}{r}\right)&\text{for}\hskip
5.69046ptR_{m}<r<R_{o}\end{dcases}$ (36)

respectively, where the dimensions are given in figure 5(a) and $Re=50$.

(a)

(b)

Figure 5: a.) Schematic of the co-annular validation study domain. b.)
Schematic of the heated channel validation study domain periodic in $x$ and
$z$ directions.

The figures 6 a,b,c, and d, show the
$L_{2}=\sum\left({\phi_{exact}-\phi}\right)^{2}/\sum{\phi_{exact}^{2}}$ and
$L_{\infty}=\text{max}|\phi_{exact}-\phi|$ norms after reaching the
convergence criteria of $L_{2}^{n+1}-L_{2}^{n}|\leq 1\cdot 10^{-9}$ ($\approx
20000$ iterations at $dt=1\cdot 10^{-4}$), and a steady state.

(a)

(b)

(c)

(d)

Figure 6: a.) The $L_{2}$ error norm for temperature and b.) the $L_{\infty}$
error norm for temperature both for $\kappa_{s}/\kappa_{f}=9$. c.) The $L_{2}$
error norm for temperature, and d.) The $L_{\infty}$ error norm for
temperature, both for $\kappa_{s}/\kappa_{f}=900$. All four plots show second
order convergence of error, $---$ is a line with LOG-LOG slope of -2 and
$\circ-\circ-\circ$ are the error norms, they are both plotted against the
number of nodes, $N=Nx=Ny$ in the simulation.

The time series show in figure 7 gives good indication that even in transience
the CHT equations are still holding the time evolution of the temperature
profiles to the sharp boundary at the discontinuity of $\kappa$. Even with
$\kappa_{s}>>\kappa_{f}$, as shown in figure 7(c) the initially discontinuity
of the time series is nearly square and the solution remains stable, even in
the points close to the interface. During the development of the CHT algorithm
defined in section 4.4 experiments were initially run purely explicitly and it
was noted that as the ratio $\kappa_{s}/\kappa_{f}$ was increased so was the
temporal stiffness of the time evolution equation making large $\kappa$
rations slow and unwieldy. The edition of the CHT boundary conditions to an
implicit solver seems to have mitigated the issues of small $dt$ for large
$\kappa_{s}/\kappa_{f}$.

For both the plots in figures 7(b) and 7(d) the final steady state lays
directly over the analytical solution given in equation 36, showing very
little deviation at any point including interfacial points. Two different were
used to find the final steady state. Initiating the temperature fields to
equation 36 and waiting until the error norms converged and letting the fields
pass through transience to the steady state solution. Both simulation types
produced indistinguishable error plots (figure 6), that were second order
convergent for $L_{\infty}$, and $L_{2}$.

(a)

(b)

(c)

(d)

Figure 7: a.) The time series for the temperature profile plotted with the
analytical solution, b.) and the final temperature profile plotted against the
analytical solution for$\kappa_{s}/\kappa_{f}=9$. c.) The time series for the
temperature profile plotted with the analytical solution, d.) and the final
temperature profile plotted against the analytical solution
for$\kappa_{s}/\kappa_{f}=900$, lines have been added to the time series to
help distinguish each time step. $-----$ is the analytical solution
$\circ-\circ-\circ$ are the temperature profiles increasing with time in
direction of the arrow. All plots shown are from the $N_{x}=256$, $N_{y}=256$
simulations with every two hundredth point removed from the plot to make it
easier to see individual points.

###  5.2 Turbulent Heated Channel Validation

The second numerical experiment for validation of the CHTIB is a turbulent
heated channel originally put forward in [8] with the addition of CHT added in
[22]. The geometry of the heated channel domain follows the simple geometry in
figure 5. The normalization in previous works, [9], [7], and [22], was not
applicable to the momentum and pressure solver used in this work so the
general equations, here, were normalized as put forward in equations 30-33.
With the exception that the normalized groups are now based on the channel
half width $L_{y}/2$ and the skin friction velocity at the wall,
$u_{\tau}=\sqrt{\mu\left(du/dy\right)_{wall}}$,

$\begin{split}\textbf{u}^{*}=\frac{\textbf{u}}{u_{\tau}}\end{split}\hskip
14.22636pt\begin{split}\textbf{x}^{*}=2\frac{\textbf{x}}{L_{y}}\end{split}\hskip
14.22636pt\begin{split}t^{*}=2t\frac{u_{\tau}}{L_{y}}\end{split}\hskip
14.22636pt\begin{split}P^{*}=\frac{P}{\rho u_{\tau}^{2}}\end{split}\hskip
14.22636pt\begin{split}Re_{u_{\tau}}=2\frac{u_{\tau}L_{y}}{\nu},\end{split}$
(37)

(a)

(b)

Figure 8: a.) The root-mean-square $\left(rms\right)$ of the fluctuating
temperature inside the heated wall b.) The $\left(rms\right)$ of the
fluctuating temperature inside the core flow. For both a.) and b.) the wall
units are those taken from [22], as well as the scaling of the temperature
field inside the wall and the core flow. Current results, $($——$)\hskip
2.84544pt$and results from [22], $\circ-\circ-\circ$.

d and the addition of a source term to equation 33,

$\partial_{t}\phi^{*}=\frac{1}{RePr}\frac{\alpha_{s}}{\alpha_{f}}\boldsymbol{\nabla}^{*}\cdot\left(\boldsymbol{\nabla}^{*}\phi^{*}\right)+\frac{\alpha_{f}}{\alpha_{s}}\frac{1}{Re_{u_{\tau}}D_{\text{wall}}}.$
(38)

where $\alpha_{s}\alpha_{f}=4$, and $Pr=7$ and the term to the to the left on
the right hand side of equation 38 is the heating source term. The boundary of
the channel domain at $y=\pm\left(\frac{L_{y}}{2}+D_{\text{wall}}\right)$ is
kept adiabatic. At $\Gamma$ the ensemble average of temperature, in time and
space, is $\overline{\phi}=0$ with the no slip condition for velocity and a
constant pressure gradient. The mesh resolution and stretching are, $Lx=64$,
$Ly=176$, $Lz=64$, and $\delta y^{+}=0.029\leftrightarrow 2.95$. The channel
is simulated until the momentum fields reach a statistical steady state where
the turbulent production is balanced by the dissipative terms [35].

in The results from [22] are scaled with the normalization used in equation 37
and compared with the results from this work in figure 8.

##  6 Conclusion

Analysis of the CHTIB algorithm shows second order accuracy with a sharp and
distinct definition between the boundaries of the fluid and solid domain. With
implicit implementation a large rage of $\kappa_{s}/\kappa_{f}$ can be
simulated without stability or accuracy concerns. Testing to larger Reynolds
numbers with 3D chaotic flow also shows that it is correlated with previous
studies giving accurate results temperature field results for large $Pr$. The
momentum IB also shows that there was little detectable effect on the
subsequent scalar field through the convective terms. The consequence is
global second order error for the scalar fields and the first order momentum
statistics

The addition of a level set field to define $\Gamma$ is a simple way to store
geometric data for the cut-cell IB geometries used by both the momentum IB and
the CHTIB, with a low computational and front end coding costs. With both the
CHTIB and momentum IB activated channel simulations for equal numbers of grid
points show a $7-10\%$ decrease in computational efficiency, overall, but
consistency with core scale up was little effected.

##  7 Acknowledgments

Our funding through the grants NSF CBET-0967857 and NASA NNX11AM07A. Professor
Olivier Desjardins for the use of his NGA flow solver, and Dr. Michael Wright
and Dr. Nagi Mansour for all their help at NASA AMES.

## References

## References

  * [1] P. Moin, K. Mahesh, Direct numerical simulation: A tool in turbulence research, Annu. Rev. Fluid Mech. 30 (1999) 39–578.
  * [2] J. Kim, P. Moin, R. Moser, Turbulent statistics in fully developed channel flow at low Reynolds number, J. Fluid Mech. 177 (1987) 133–166.
  * [3] R. Moser, J. Kim, N. Mansour, Direct numerical simulation of turbulent channel flow up to r=559, Physics of Fluids 11 (4) (1999) 943–945.
  * [4] Y. M. Khabakhpasheva, Experimental Investigation of Turbulent Momentum and Heat Transfer in the Proximity of the Wall, Heat Transfer 1986, Proceedings of The 8th International Heat Transfer Conference.
  * [5] M. Teitel, R. Antonia, Heat transfer in fully developed turbulent channel flow: comparison between experiment and direct numerical simulations, International Journal of Heat and Mass Transfer 36 (6) (1993) 1701–1706.
  * [6] Y. L. Sinai, A Wall Function for the Temperature Variance in Turbulent Flow Adjacent to a Diabatic Wall, J. Heat Transfer 109 (1987) 861–865.
  * [7] H. Kawamura, K. Ohsaka, H. Abe, K. Yamamoto, DNS of turbulent heat transfer in channel flow with low to medium-high Prandtl number fluid, Heat and Fluid Flow 19 (1998) 849–857.
  * [8] N. Kasagi, A. Kuroda, M. Hirata, Numerical Investigation of Near- Wall Turbulent Heat Transfer Taking Into Account the Unsteady Heat Conduction in the Solid Wall, Journal of Heat Transfer 111 (1989) 385–392.
  * [9] N. Kasagi, Y. Tomita, A. Kuroda, Direct Numerical Simulation of Passive Scalar Field in a Turbulent Channel Flow, Journal of Heat Transfer 114 (1992) 598–606.
  * [10] J. Kim, P. Moin, Transport of Passive Scalars in a Turbulent Channel Flow, Turbulent Shear Flows 6 6 (1989) 85–96.
  * [11] S. Kang, G. Iaccarino, F. Ham, DNS of buoyancy-dominated turbulent flows on a bluff body using the immersed boundary method, J. Comp. Phys. 228 (9) (2009) 3189–3208.
  * [12] S. Kang, An Improved Immersed Boundary Method for Computation of Turbulent Flows with Heat Transfer, Ph.D. thesis, Stanford University (2008).
  * [13] K. Nagendra, D. Taftia, K. Viswanathb, A new approach for conjugate heat transfer problems using immersed boundary method for curvilinear grid based solvers, J. Comp. Phys. 267 (2014) 225––246.
  * [14] C. Peskin, Flow patterns around heart valves: A digital computer method for solving the equations of motion, PhD Thesis, Albert Einstein College of Medicine.
  * [15] J. Mohd-Yusof, Combined immersed boundaries/B-splines methods for simulations of flows in complex geometries, Center for Turbulence Research Annual Research Briefs NASA Ames/Stanford University.
  * [16] J. Jongho Lee, Kim, H. Choi, K. Yang, Sources of spurious force oscillations from an immersed boundary method for moving-body problems, Journal of Computational Physics 230.
  * [17] T. Ye, R. Mittal, U. HS, W. Shyy, An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries, J. Comp. Phys. 156 (1) (1999) 209–240.
  * [18] R. Mittal, H. Dong, M. Bozkurttas, F. Najjar, A. Vargas, A. Von Loebbecke, A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries, J. Comp. Phys. 227 (10) (2008) 4825–4852.
  * [19] M. Kirkpatrick, S. Armfield, J. Kent, A representation of curved boundaries for the solution of the navier–stokes equations on a staggered three- dimensional cartesian grid, J. Comput. Phys. 184 (2003) 1–36.
  * [20] M. Meyer, A. Devesa, S. Hickel, X. Hu, N. Adams, A conservative immersed interface method for large-eddy simulation of incompressible flows, Journal of Computational Physics 229 (2010) 6300–6317.
  * [21] O. Desjardins, G. Blanquart, G. Balarac, H. Pitsch, High order conservative finite difference scheme for variable density low Mach number turbulent flows, J. Comp. Phys. 227 (2008) 7125–7159.
  * [22] I. Tiselj, R. Bergant, B. Mavko, I. Bajsic, G. Hetsroni, DNS of Turbulent Heat Transfer in Channel Flow With Heat Conduction in the Solid Wall, Journal of Heat Transfer 123 (2001) 849–857.
  * [23] S. Osher, J. Sethian, Fronts Propagating with Curvature Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations, J. Comput. Phys. 79 (1988) 12–49.
  * [24] O. Desjardins, J. McCaslin, M. Owkes, P. Brad, Direct numerical and large-eddy simulation of primary atomization in complex geometries, Atomization and Sprays 23 (2013) 1001–1048.
  * [25] Y. Morinishi, T. S. Lund, O. V. Vasilyev, P. Moin, Fully Conservative Higher Order Finite Difference Schemes for Incompressible Flow, J. Comp. Phys. 143 (1998) 90–124.
  * [26] Y. Morinishi, O. V. Vasilyev, T. Ogi, Fully conservative finite difference scheme in cylindrical coordinates for incompressible flow simulations, J. Comp. Phys. 143 (2004) 686–710.
  * [27] S. Osher, F. Solomon, Upwind difference schemes for hyperbolic systems of conservation laws, Mathematics of Computation 158 (1982) 339–374.
  * [28] K. Henrick, T. Aslam, J. Powers, Mapped weighted essentially non-oscillatory schemes: Achieving optimal order near critical points, J. Comp. Phys. 207 (2005) 542–567.
  * [29] J. Kim, P. Moin, Application of a fractional-step method to incompressible Navier-Stokes equations, J. Comp. Phys. 59 (1985) 308–323.
  * [30] P. M. H. Choi, Effects of the computational time step on numerical solutions of turbulent flow, J. Comp. Phys. 113 (1) (1994) 1–4.
  * [31] C. Pierce, Progress-variable approach for large eddy simulation of turbulent combustion, Ph.D. thesis, Flow Physics and Computation Division, Dept. Mech. Eng., Stanford Univ (2001).
  * [32] Z. L. Chen, S. Hickel, A. Devesa, J. Berland, N. Adams, Wall modeling for implicit large-eddy simulation and immersed-interface methods, Theoretical and Computational Fluid Dynamics 28 (2014) 1–21.
  * [33] T. Ye, R. Mittal, H. Udaykumar, W. Shyy, An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries, Journal of Computational Physics 156 (1999) 209–240.
  * [34] X. Hu, B. Khoo, N. Adams, F. Huang, A conservative interface method for compressible flows, J. Comput. Phys 219 (2006) 553–578.
  * [35] S. B. Pope, Turbulent Flows, Cambridge University Press, 2000.