Bruno Bierbaumer
email: bierbaumer@sec.in.tum.de

Julian Kirsch

Thomas Kittel

Aurélien Francillon

Apostolis Zarras

Smashing the Stack Protector for Fun and Profit

Software exploitation has been proven to be a lucrative business for cybercriminals. Unfortunately, protecting software against attacks is a long-lasting endeavor that is still under active research. However, certain software-hardening schemes are already incorporated into current compilers and are actively used to make software exploitation a complicated procedure for the adversaries. Stack canaries are such a protection mechanism. Stack canaries aim to prevent control flow hijack by detecting corruption of a specific value on the program's stack. Careful design and implementation of this conceptually straightforward mechanism is crucial to defeat stack-based control flow detours. In this paper, we examine 17 different stack canary implementations across multiple versions of the most popular Operating Systems running on various architectures. We systematically compare critical implementation details and introduce one new generic attack vector which allows bypassing stack canaries on current Linux systems running up-to-date multi-threaded software altogether. We release an open-source framework (CookieCrumbler ) that identifies the characteristics of stack canaries on any platform it is compiled on and we propose mitigation techniques against stack-based attacks. Although stack canaries may appear obsolete, we show that when they are used correctly, they can prevent intrusions which even the more sophisticated solutions may potentially fail to block.

Introduction

Buffer overflow vulnerabilities are as old as the Internet itself. In 1988, the Morris Worm was one of the first malware discovered in public that leveraged this vulnerability [START_REF] Streak | The Morris Worm: A Fifteen-Year Perspective[END_REF]. Since then, many security breaches can be linked to successful exploitation of buffer overflows, which denotes that the problem is far from being solved. As a matter of fact, the Mitre Corporation lists more than eight thousand Common Vulnerabilities and Exposures (CVE) entries that contain the keyword "buffer overflow" 4 . A significant portion of these vulnerabilities is comprised by the so-called stack-based buffer overflow bugs [START_REF] One | Smashing the Stack for Fun and Profit[END_REF]. This is due to the application's stack inherent property of mixing user-controlled program data together with control flow related data; thus allowing an attacker to overwrite control flow related parts of the stack.

Due to this well-known weakness, Cowan et al. [START_REF] Cowan | StackGuard: Automatic Adaptive Detection and Prevention of Buffer-Overflow Attacks[END_REF] propose a technique named Stack Smashing Protection (SSP). The idea behind SSP is to detect stack-based control flow hijacking attempts by introducing random values (so-called canaries) to the stack that serve as a barrier between attacker-controlled data and control flow relevant structures. After a function finishes executing, a canary-named after the coal miner's canaries which were used to detect presence of gas-is checked against a known "good" value stored in a safe location. Only if the canary maintains its original value, execution continues. This mitigation technique has been present in compilers for more than 10 years and is now a counter measure supported by major compilers [START_REF]Using the GNU Compiler Collection[END_REF][START_REF]Clang's Documentation[END_REF][START_REF]Compiling a C/C++ Program[END_REF].

Recently, more advanced techniques have been proposed to prevent buffer overflow attacks including Code Pointer Integrity (CPI) [START_REF] Cowan | Buffer Overflows: Attacks and Defenses for the Vulnerability of the Decade[END_REF][START_REF] Kuznetsov | Code-Pointer Integrity[END_REF] and Control Flow Integrity (CFI) [START_REF] Abadi | Conference on Computer and Communications Security (CCS)[END_REF]. Both ideas encircle the concept of protecting the control flow from being hijacked. These advanced techniques have created the illusion that stack canaries are nowadays obsolete. However, both techniques consider non-control-flow diverting attacks to be out of scope. As we discuss later, this is an underestimated attack that can be successfully countered by stack canaries [START_REF] Chen | Non-Control-Data Attacks Are Realistic Threats[END_REF]. While introduced almost twenty years ago, stack canaries are still one of the most widely deployed defense mechanisms to date [START_REF] Szekeres | Eternal War in Memory[END_REF] and are, as we will show, a necessary complement to other more recent modern buffer overflow mitigation mechanisms. As a matter of fact, all modern compilers support stack canaries.

In this paper, we show that stack canaries, even in combination with more advanced techniques, are not a silver bullet. We find that due to inconsiderate implementation decisions, stack canaries themselves are vulnerable to buffer overflow attacks; ironically the same type of attack that they are supposed to protect against. To demonstrate this, we first implement a framework (CookieCrumbler ) which is able to identify the characteristics of stack canaries on various modern operating systems, independently from the CPU architecture used. We then run CookieCrumbler on 17 different combinations of Operating Systems (OSes), C standard libraries, and hardware architectures. We run CookieCrumbler against both seemingly old, but still widely used and supported, OSes as well as the most recent versions. The extracted CookieCrumbler results enable us to introduce a new attack based on the observation that the canary reference values are not always stored at a safe location. This allows an attacker to overwrite the control-flow relevant data structures and the current reference value at the same time and thereby reliably bypass SSP.

In summary, we make the following main contributions:

-We propose CookieCrumbler , a framework to automate the identification of the characteristics of SSP implementations. -We evaluate CookieCrumbler against state of the art operating systems and libraries, and discover weaknesses in multiple SSP implementations. -We introduce a novel attack vector to exploit those vulnerabilities.

-We propose mitigations techniques to harden SSP implementations.

Background and Related Work

Stack Smashing Protection

The idea to guard certain parts of the executable's stack dates back to 1998 [START_REF] Cowan | StackGuard: Automatic Adaptive Detection and Prevention of Buffer-Overflow Attacks[END_REF][START_REF] Cowan | Protecting Systems From Stack Smashing Attacks With StackGuard[END_REF]. The concept is to protect control flow related information on the stack using a so-called stack canary or stack cookie: a random value placed between the user-controllable data and the return pointers on the stack during stack setup phase in the function prologue. The mechanism is implemented synchronously with the control flow: after function execution, once the control flow returns to the caller, the cookie value is checked against a known "good" value. Only if there is a match between the two values, the stack frame is cleaned up and the control flow is allowed to return to the caller. Several years were needed for StackGuard to be integrated in the mainline GCC distribution [START_REF] Wagle | Stackguard: Simple Stack Smash Protection for Gcc[END_REF].

Attackers may try to evade StackGuard by embedding the canary in the data used during the overflow (i.e., canary forgery). Cowan et al. [START_REF] Cowan | Buffer Overflows: Attacks and Defenses for the Vulnerability of the Decade[END_REF] propose two methods to prevent such a forgery: terminator and random canaries. In 32-bit operating systems, a terminator canary is usually constructed using of the char representation of NULL, CR, LF, and EOF (0x000d0aff). This is because overflows often exploit unsafe string manipulation functions: as the terminator canary includes characters which are used to terminate strings, it is impossible to directly include them in a string without terminating the string operation. However, not all buffer overflows are due to unsafe string manipulation operations (e.g., read()) and the fixed terminator canary does not provide any protection in those cases. On the other hand, random canaries cannot be guessed by an attacker and is therefore the most generic approach.

Marco-Gisbert and Ripoll extend the original StackGuard concept by proposing a renewal of the secret stack canary during the fork and clone system calls [START_REF] Marco-Gisbert | Preventing Brute Force Attacks Against Stack Canary Protection on Networking Servers[END_REF]. This way, an external attacker is not able to brute-force the stack canary in scenarios where the request handling routine is forked from a server application for each request, as is typically the case for network facing applications. As an alternative, Kuznetsov et al. [START_REF] Kuznetsov | Code-Pointer Integrity[END_REF] propose secure code pointers by storing them in a safe memory region. In their work, they assume that the location of the safe region can be hidden.

In essence, in order to be effective, StackGuard relies on the following assumptions:

The cookie value placed on the stack (Can) must be unknown to the attacker. The known good value (Ref) is placed at a location in memory that is distinct from the location of Can and ideally mapped read-only. If a stack cookie value (Can) is corrupted, the program execution terminates immediately without accessing any attacker controlled data. The overflow is contiguous, starting from a buffer on the stack, and therefore does not permit an attacker to skip certain bytes in memory.

The main focus of this paper is to falsify Assumptions , , and using contiguous overflows (i.e., adhering to Assumption ).

Function Pointer Protection

Function pointers stored in writable memory at static addresses are common targets to gain control of a vulnerable program's execution. To defend against this threat, [START_REF] Cowan | Buffer Overflows: Attacks and Defenses for the Vulnerability of the Decade[END_REF] introduces code pointer protection. PointGuard [START_REF] Cowan | Pointguard: Protecting Pointers From Buffer Overflow Vulnerabilities[END_REF] is the first mechanism capable of encrypting code pointers in memory. For each process, PointGuard generates a random key during process creation. Each pointer in memory is then scrambled in memory by performing a bijective operation on the pointer using the process's specific random key.

Glibc implements this protection mechanism since 2005 [START_REF] Drepper | Pointer Encryption[END_REF] by using the PTR _ [DE]MANGLE macros. On the other hand, the Windows run-time, provides similar functionality with the Rtl[En|De]codePointer API call since XP SP2 (2004). Both implementations use very similar algorithms to encipher pointers: a logical bit rotation combined with an xor (⊕) involving the per-process random secret (rand). For instance, the 64-bit Windows run-time implements the following two equations for pointer protection:

ptr enc = ror64(ptr orig ⊕ rand, rand) (1)
ptr orig = rol64(ptr enc , rand) ⊕ rand (2)
On Linux (with glibc), the situation is very similar, except that a constant is used as the number of digits to rotate (0x11 is actually 2 • sizeof(void * ) + 1):

ptr enc = ror64(ptr orig ⊕ rand, 0x11) (3)
ptr orig = rol64(ptr enc , 0x11) ⊕ rand (4)
Most notably, the main difference in the two implementations is that on Windows the Rtl[En|De]CodePointer retrieves the value rand from the kernel, whereas glibc on Linux stores the pointer guard in user space in the Thread Control Block (TCB). Cowan et al. [START_REF] Cowan | Pointguard: Protecting Pointers From Buffer Overflow Vulnerabilities[END_REF] state that the PointGuard key has to be stored on its own page once it is initialized to protect the key against information leakage. As can be seen, this assumption is not met by all implementations.

Attacks Against Stack Canaries

An adversary may attempt to attack the stack canary mechanism itself in order to successfully exploit a program. Strackx et al. [START_REF] Strackx | European Workshop on System Security (EUROSEC)[END_REF] analyze the security promises made by randomization based buffer overflow mitigation systems, such as the ones described above. They conclude that a vulnerable program offering both a buffer overread and a buffer overflow can be easily attacked. However, their work misses the experimental evaluation of the success rate of such an attack. Ding et al. [START_REF] Ding | Android Low Entropy Demystified[END_REF] reveal weaknesses in the StackGuard implementation used in Android 4.0: the source of randomness used for the stack canaries is only initialized once at OS boot and then used for every application on the system. In addition, the created canary is predictable as the state used to initialize the canary only depends on randomness available at kernel boot-up.

Dynamic Canary Randomization [START_REF] Hawkins | Dynamic Canary Randomization for Improved Software Security[END_REF] attempts to defend attacks targeting stack canaries. This technique re-randomizes all active stack canaries during run-time so the attackers cannot reuse the knowledge they gained while leaking memory from an earlier execution of the attacked process. While this approach might help against attacks that read the canary and then use the gained knowledge in a separate step, it is ineffective against the attack introduced in this work.

Thread Control Block

Modern OSes contain a dedicated data structure, called TCB, to store information about the environment that the current thread is executing in. The data stored in the TCB varies depending on OSes and thread library implementations. For instance, on Windows this data structure is named ThreadInformationtBlock and contains information about a thread's Structured Exception Handling (SEH) chain, its associated Process Control Block (PCB), and a pointer to Thread Local Storage (TLS). The TCB is accessed either using a library function or a designated register that improves speed. For example, glibc on Linux x86_64 uses the fs register as the base address of the TCB. Intel provides Model Specific Register s (MSRs) to override fs and gs segment base addresses, effectively enabling 64-bit OSes to access the TCB in a fast way. This is achieved by prefixing any load or store operation with the fs segment register.

Both SSP (with StackGuard) and Function Pointer Protection (with Point-Guard) belong to the standard set of defense mechanisms and are highly adopted in practice. However, both mechanisms require the storage of their respective random reference keys (Ref). This is where the TCB becomes relevant, in context of our work: In some versions of the compiler/standard library, it is the TCB that contains the reference keys for both mechanisms. Therefore, both mechanisms can be attacked if the data contained in the TCB can be overwritten, as shown in this paper.

Modern Defense Mechanisms

In this study, we explicitly concentrate on concrete implementations rather than theoretic contributions. We therefore focus on defense mechanisms that are (i) available in current (2018) compilers, (ii) production ready, and (iii) deployed in current operating systems. This leaves us with a very narrow set of mechanisms. In fact, we rarely find academic publications implementations that reach a mature state PointGuard [START_REF] Cowan | Pointguard: Protecting Pointers From Buffer Overflow Vulnerabilities[END_REF], and StackGuard [START_REF] Cowan | Protecting Systems From Stack Smashing Attacks With StackGuard[END_REF].

There is a trend on maintaining the integrity of an application's control flow at runtime. CFI is achieved by ensuring the integrity of forward and backward edges in the control flow graph. As we focus on stack based exploitation techniques targeting control flow information related to the backward edge, we only consider the backward edge validation relevant. Backward edge validation is typically done using a shadow stack [START_REF] Abadi | Conference on Computer and Communications Security (CCS)[END_REF]. One production ready implementation is SafeStack [START_REF] Kuznetsov | Code-Pointer Integrity[END_REF], which we inspect in Section 6 to understand its relationship to stack canaries. Algorithm 1: Algorithm used to measure empirical features.

Dissecting Implementation Choices

In this section we define five qualitative and five empiric features which we use to systematically evaluate choices made in canaries implementation.

Qualitative Features

We identify key features of stack canaries by studying the source code of their implementations-if available-or reverse engineering the functionality in their binary format. Assumption (immediate termination) is another claim that can only be verified in a qualitative manner. To find the quantity of code executed after the canary corruption is detected, we introduce the notion of noisiness of the failure handler. To estimate the Noise level, we count function invocations that are triggered from the point where execution enters the cookie verification failure handler until the point where the application terminated. We also manually check the number of variables that are read from the corrupted memory region (e.g. the stack), and whether the handler executes in user or kernel mode, which we denote by Current Privilege Level (CPL).

Empirical Features

To reason about potential attack targets, we retrieve basic information about the application's memory layout. For each OS and C library pair, we run a test program which follows Algorithm 1. The program measures the distance (in terms of their addresses) between each user-controllable types of memory. This distance measurement is an important information. Indeed, the closer the Ref value is from a user controllable memory the easier it will be for an attacker to overwrite this reference value, and therefore to be able to corrupt the canary without being detected.

More precisely we measure spatial distances (∆) between the reference value (Ref) and:

1. ∆ loc : a variable allocated on the stack of the function. We then compute the range of non-contiguous bytes in ∆ x . If this range is not mapped as a contiguous writable memory, an overflow from this variable will trigger a page fault before reaching Ref.

5. W(∆ x ): number of contiguously mapped writable bytes in ∆ x .

CookieCrumbler

We implemented the CookieCrumbler framework to evaluate those features.

From a high-level perspective, CookieCrumbler is a direct implementation of Algorithm 1 in C. When compiled and executed on a system, CookieCrumbler will thoroughly analyze the implementation of stack canaries. For this purpose, semantic knowledge about the exact location of Ref has to be added to the program. For instance, on x86_64, Ref is located within the TCB at offset 0x28.

We include this information for all the environments presented in Section 4.

The core of Algorithm 1 is to retrieve the deltas ∆ loc , ∆ glo , ∆ dyn , and ∆ tls . To obtain the respective reference point in memory, we use (i) a stack local variable, (ii) a variable with the static keyword, (iii) the pointer value returned by malloc, and (iv) a variable with the __ thread keyword (on UNIX) or __ declspec(thread) (on Windows). Threads are created by calls to the functions pthread _ create (on UNIX) or CreateThread (on Windows). To determine W(∆ x ), we use signal handling on UNIX (catching SIGSEGV on a contiguous byte-by-byte write) and the function IsBadWritePtr on Windows.

After successful execution, CookieCrumbler generates a set of memory locations, deltas, and number of writable bytes for the main-and the sub-threads of a threaded application, respectively. A thorough analysis of these results can reveal potential vulnerabilities in the implementation of stack canaries. The source code and the measured data can be found online. 5

Smashing the Stack Protector

We run CookieCrumbler on various OSes with different C standard libraries. Apart from up-to-date version of the C runtime libraries, we also run CookieCrumbler on older libc versions that are currently still distributed in the stable branches of commonly used Linux distributions. For more details refer to Table 1.

Qualitative Results

Surprisingly, the qualitative features we examined look very homogenous. We therefore first explain the most common observations and then discuss special cases. Unexpectedly, we found that almost none of the tested implementations changes Can across different function invocations within the context of one given thread. The only exception to this rule constitutes the Windows family of operating systems, for which Can is chosen as Ref ⊕ rbp when the rbp register is used as stack frame pointer and Can = Ref ⊕ rsp otherwise.

As indicated by the literature [START_REF] Marco-Gisbert | Preventing Brute Force Attacks Against Stack Canary Protection on Networking Servers[END_REF] we also observed Ref (and consequently Can for all stack frames) to remain static across fork invocations on all UNIX operating systems. On this particular point, comparison with Windows is impossible as the fork system call is not supported by Windows operating systems family.

In nearly all cases the failure handler executes in user space (the privilege level CPL is 3). The only exceptions to this rule are Windows 8 and newer, which implement the special interrupt number 0x29 ( _ KiRaiseSecurityCheckFailure) for this purpose. When this interrupt handler is called, the program is terminated without accessing any of the potentially corrupted memory in userspace. Windows can fall back to the old user-space failure routine if a call to IsProcessorFeaturePresent(PF _ FASTFAIL _ AVAILABLE) returns zero.

On Windows OS versions newer than 7, the Noise level is the lowest, as they support an interrupt specifically designed for this purpose. Older versions call 8 functions in kernel32.dll and collect information about the current register state before terminating (TerminateProcess) the application with return code 0xc0000409 (Security check failure or stack buffer overrun). OpenBSD, when detecting a corrupt stack canary, infers the program's name from a (safe) location in the global variable section of the currently loaded standard library and prints one line of information into the system log. Linux's C standard libraries implement __ stack _ chk _ fail in different ways: musl libc does not provide any output and terminates execution using a hlt instruction, accounting for a minimal Noise level. diet libc prints a static error message and terminates the program with an exit syscall. Bionic logs a static message, which required to allocate dynamic memory, and finally terminates the program via a SIGABRT. The Noise level culminates on Linux with glibc prior to version 2.26: where we measured that the __ stack _ chk _ fail function performs as many as 69 calls to other functions, dispatching at least three calls using (PointGuard protected) writable global static function pointers to create a stack trace by unwinding the attacker controlled stack before exiting the process. More importantly, glibc prints the program name fetched from the argv array on the stack, which is a potentially attacker-controlled location creating an arbitrary memory leak primitive. This behavior (assigned CVE-2010-3192) was finally fixed in glibc version 2.26 in August 2017.

Empirical Results

We classify our data points into three categories:

1. The vulnerable implementations satisfying ∆ loc > 0 and W (∆ loc ) = 100.0% are marked in . Here, a long buffer overflow on the stack allows for a complete stack canary bypass as Can and Ref can be overwritten at the same time. 2. The weak implementations satisfying W (∆ ) = 100.0% with ∆ = ∆ loc are marked in . This requires an attacker to not only overflow a data structure located in the memory segment next to Ref (maybe even in reverse direction), but also to get control of the execution flow by overwriting a buffer on the stack before the function containing the first vulnerability returns. 3. The secure implementations satisfying W (∆) = 100.0% are marked as in Table 1. These implementations do not offer the possibility to overwrite Ref in memory and therefore are secure against the attack presented in this work.

In essence, Categories 1 and 2 violate the Assumption .

Introduced Attack Vectors

We now discuss the practical implications for application security. For clarity we omit the discussion of weak implementations, as an attacker would always need more than one buffer-overflow in a vulnerable application to gain advantage of the situation. For this, we assume an adversary who is capable of triggering a buffer-overflow of suitable size on the stack. As such, we discuss possible attack vectors in two different scenarios depending on the threading model the target executable uses.

Forking. In a forking environment, the whole address space of the target binary is duplicated, including all Can and Ref values contained in memory. When an attacker is able to obtain information about one of the forked processes this renders randomness based countermeasures ineffective as all forked applications share the same randomness: ASLR becomes predictable [START_REF] Shacham | On the effectiveness of address-space randomization[END_REF] as well as all cookie values. Assuming an attacker is allowed to restart communication with the vulnerable application an oracle can be created as follows: the attacker overwrites a stack canary byte by byte and observes whether the application at the other end crashes. Only one out of 2 8 possible byte values will allow the application to continue execution. This effectively increases the chance of guessing the stack canary from (2 8 ) 8 = 2 64 to (2 8 ) • 8 = 2 11 in the worst case-implying a more than significant difference in both attack duration as well as probability of success. This attack vector has already been discussed by researchers [START_REF] Marco-Gisbert | Preventing Brute Force Attacks Against Stack Canary Protection on Networking Servers[END_REF][START_REF] Hawkins | Dynamic Canary Randomization for Improved Software Security[END_REF][START_REF] Bittau | Hacking blind[END_REF]. Note that a similar technique can be used to infer certain pointer values residing in the attacked application's stack frame.

Threading. On multi-threaded applications the insights from CookieCrumbler can be used in two ways. (1) If the attacker can write null bytes and the application is mapped at a static address in memory: all vulnerable implementations stack canaries can be completely bypassed by overwriting Can and Ref with the same value chosen by the attacker. As all program addresses are known, this case directly reduces to an ordinary Return Oriented Programming (ROP) attack.

(2) If the attacker is not allowed to write null bytes or the application's code section is not mapped at a static address (e.g., Position Independent Executable (PIE)): the attack can still succeed on Linux with glibc. The attacker will target the PointGuard value, which is also stored in the TCB, directly following Ref.

Equations 3 and 4 show that in PointGuard any protected pointer is first rotated by a fixed number of digits and then xored with the PointGuard value (i.e., an attacker controlled number) in the considered setting.

The function in charge of terminating the program after a failed stack cookie check in glibc eventually ends up demangling a pointer to pthread _ once. It is obvious that by the simple arithmetic used during pointer demangling, the attacker can detour the execution flow by a fixed offset to this function. From here on, no generic attack vector exists, but we want to point out that there are code paths in glibc that execute the assembly-equivalent of execve("/bin/sh") [START_REF]One Gadget[END_REF], which constitute valuable attack targets in our case. The likelihood of this attack succeeding heavily depends on the memory layout imposed by the dynamic loader on libraries. In our experiments we never observed a distance greater than 2 24 between pthread _ once and a gadget that eventually lead to remote code execution.

Impact

The tested Linux-based platforms (Android, Arch Linux, Debian, and Ubuntu) can be clustered into two different categories. Architectures which have dedicated TLS access registers (x86, x86_64, s390x, and PowerPC) that store the Ref in the TCB and architectures without a direct register access to the TLS (ARM). We have also analyzed the source code of glibc and categorized further architectures as TLS-based stack canary implementations: IA64, SPARC, and TILE. While we expect that our results can be extended to those architectures we did not had access to such hardware and did not include them in Table 1.

The TLS-based SSP for all the libc implementations are vulnerable to our attack in a multi-threaded environment by overwriting the Ref via a stack-based buffer overflow 6 . SSP implementations where the Ref is located in the Global section are more robust as it cannot be modified by an buffers overflow on the stack. This result can be seen in the loc column in Table 1. Our evaluation also shows that most implementations fail to separate other data regions from the location of Ref. This might be exploitable if the program uses thread-local variables. If one of the variables can be overflown, an attacker may overwrite the reference canary Ref. In this case, the attacker needs two overflows (to change both Ref and Can). This is a difficult attack which also affects single-threaded applications, and therefore less critical issue.

Interestingly, diet libc defaults to storing the reference canary in the TLS, even if the application is not multi-threaded. Thus also the main thread stack is adjacent to the used TLS. Note that the main thread's stack and its TLS region are separated in the other implementations. This effectively breaks SSP for diet libc. Also, we point out that SSP can be bypassed for multi-threaded applications in all libc implementations.

Windows, macOS, and BSD derivatives store the reference cookie in the .bss section. Hence, they are not vulnerable to our overflow attack. However, column glo in Table 1 shows that storing the reference stack cookie in the .bss region might open up a vulnerability. On Windows and FreeBSD, the stack canary is located in front of the global variables. Thus, the value might get overwritten by an overflow running towards lower addresses, which is less common yet not impossible. Only macOS, OpenBSD, and Android (on architectures without TLS based cookies, e.g., ARM) succeed to separate the reference cookie from all other memory regions. As OpenBSD, at the time of writing, is lacking a compiler with support for thread-local variables, it is not conducted in our experiments.

To get an overview how realistic the described attack is, we analyze the binaries installed on a vanilla Debian Jessie installation. About 40% of those programs depend on pthreads, which leaves them potentially vulnerable to our attack. Server applications, like web servers, often rely on threading to handle multiple clients at the same time and are particularly subject to such attacks.

Attack Mitigations

Re-randomizing Ref on process creation (e.g., after forking) is a promising idea to increase canary entropy, as demonstrated by RenewSSP [START_REF] Marco-Gisbert | Preventing Brute Force Attacks Against Stack Canary Protection on Networking Servers[END_REF]. This approach mitigates our attack partially, but we also propose to modify the thread library to randomize Ref for each thread.

Frantzen et al. [START_REF] Frantzen | StackGhost: Hardware Facilitated Stack Protection[END_REF] argue to relocate Ref to the PCB data structure, but unfortunately this introduced more deficiencies. We extend this idea by proposing the generation of per-function stack cookie values by xoring the static canary with the current stack pointer value to borrow randomness from mmap. Similarly, we can xor Ref with the return address of the protected function. However, this mitigation is only effective for scenarios where the code segment of the protected function is mapped at randomized addresses. 7Handlers running in a corrupted program context should strive to quit execution as fast as possible. Glibc's __ stack _ chk _ fail handler is a bad counter-example. It passes control through several layers of code which uses attacker controlled values from the stack. This opens the possibility for further exploitation. Clearly, the approaches taken by Microsoft Visual C (MSVC) and musl libc are preferablethe handler quits as fast as possible and, in case of MSVC, any reasoning about the crashed program's state (if at all) is performed using run-time data from the OS's kernel only.

Finally, the TCB must not be mapped adjacently to any memory structure that contains user-controllable buffers. The most direct way to achieve this is the introduction of a mandatory guard page mapped with no access protection at the bottom of the stack. Note that even though glibc's pthread implementation apparently offers such functionality (pthread _ attr _ setguardsize), it is not automatically turned on by software intending to use threads and even more importantly only offers a mechanism to map a guard page on the top of the stack.

Improving Sophisticated Protection Mechanisms

To highlight how stack canaries can improve application security, we consider the C program in Figure 1. Depending on the mitigation mechanisms added when compiling this program, the authentication bypass can be trivially triggered. We consider this example in the context of two software protection mechanisms: SafeStack: SafeStack [START_REF] Kuznetsov | Code-Pointer Integrity[END_REF] is a State-of-the-Art CPI implementation that logically separates the architectural stack into a safe and unsafe region. The safe region contains all control-flow related data while the unsafe stack contains user-controlled data (e.g., arrays). Stack Canaries: We use the standard implementation employed by LLVM. When compiling with SafeStack enabled, the variables password and admin _ hash are allocated in the unsafe stack, whereas the return addresses are in the safe stack. The stack-based buffer overflow in the auth function make bypassing the security check trivial: an attacker first overflows the password buffer and then overwrites the admin _ hash with the hash matching the provided password. When stack canaries are enabled the attack is no longer trivial. After filling the password buffer, an attacker has to overwrite Can to reach the admin _ hash. Once the auth function returns, the canary corruption will be detected and the program will terminate. The same security properties (protecting buffers of adjacent stack frames) are achieved regardless of the usage of SafeStack. To reach the admin _ hash buffer, an attacker has to overwrite the Can value. As corrupting Can should result in program termination-the fact that the return address ret is also reachable by the overflow becomes irrelevant. While SafeStack threat model does not include the corruption of non-control-flow-related data structures we argue that stack canaries can improve resistance of CPI against non-control-flow targeting attacks.

Conclusion

In this work we presented CookieCrumbler , a multi-platform framework to systematically study stack canary implementations. We discovered scenarios which are prone to a novel attack that allows bypassing State-of-the-Art stack protection mechanisms in threaded environments. In addition, we introduced new ideas for a more advanced attack that abuses the way exception routines and pointer mangling mechanisms work together. Finally, we believe this work provides systematic insight into the qualitative implementation details of stack canaries used by modern OSes and can serve as a basis for future explorations of security critical parts of the OSes and C standard libraries in use today.

  s u b m a in s u b m a in s u b m a in s u

  int auth(char * valid) { char password[32]; gets(password); return strcmp(valid, crypt(password, valid)) == 0; } int main(void) { char admin _ hash[] = "$1$01234567$b5lh2mHyD2PdJjFfALlEz1"; if (auth(admin _ hash)) puts("Welcome to the Admin Area"); }

Fig. 1 .

 1 Fig. 1. Different stack layouts of a C program exposing an authentication bypass vulnerability when using (a) SafeStack (b) SafeStack with Canaries, (c) Canaries only.

  2. ∆ tls : a variable allocated in Thread Local Storage (TLS). 3. ∆ glo : a global variable allocated in statically allocated memory. 4. ∆ dyn : a variable allocated in dynamically allocated memory.

Table 1 .

 1 Summary of problems found with CookieCrumbler .

https://cve.mitre.org/cve/search _ cve _ list.html

https://bierbaumer.net/security/cookie/

During the paper's review this issue was independently discovered by Ilya Smith: https://github.com/blackzert/aslur/

OpenBSD very recently added "RETGUARD" which is similar to our proposition https://marc.info/?l=openbsd-cvs&m=152824407931917&w=2.

Acknowledgments

The research was supported by the German Federal Ministry of Education and Research under grant 16KIS0327 (IUNO) as well as the SeCiF project within the French-German Academy for the Industry of the future


Suelen Gasparin

Julien Berger

Denys Dutykh

Nathan Mendes

An innovative method to determine optimum insulation thickness based on non-uniform adaptive moving grid

Keywords: Optimum insulation thickness, thermal insulation, numerical simulation, moving grid method, redistribution schemes, adaptive numerical methods MSC: [2010] 35R30 (primary), 35K05, 80A20, 65M32 (secondary) PACS: [2010] 44.05.+e (primary), 44.10.+i, 02.60.Cb, 02.70.Bf (secondary)

It is well known that thermal insulation is a leading strategy for reducing energy consumption associated to heating or cooling processes in buildings. Nevertheless, building insulation can generate high expenditures so that the selection of an optimum insulation thickness requires a detailed energy simulation as well as an economic analysis. In this way, the present study proposes an innovative non-uniform adaptive method to determine the optimal insulation thickness of external walls. First, the method is compared with a reference solution to properly understand the features of the method, which can provide high accuracy with less spatial nodes. Then, the adaptive method is used to simulate the transient heat conduction through the building envelope of buildings located in Brazil, where there is a large potential of energy reduction. Simulations have been efficiently carried out for different wall and roof configurations, showing that the innovative method efficiently provides a gain of 25% on the computer run time.

Introduction

In Brazil, the electricity consumption in residential, commercial and public sectors accounted for 42.8% of the Brazilian electricity consumption in 2016 [START_REF] Filho | Brazilian Energy Balance[END_REF]. The consumption of electrical energy in buildings is related to gains or losses of heat through the building envelope that is also associated with the internal loads. In particular, the use of heating, ventilation, and air conditioning (HVAC) systems can represent 40% of the total electrical energy consumption [START_REF] Abrava | [END_REF].

Due to current environmental issues, there is a great potential to reduce heating and cooling loads particularly through improvements of thermal performance of the envelope. One of the most efficient ways to reduce the transmission loads and, consequently, reduce the energy consumption is the use of an appropriated thermal insulation for the building envelope. However, in Brazil, the insulation costs are relatively high, due especially to the low demand and the lack of a compulsory regulation [START_REF] Morishita | Issues about moisture in residential buildings of Brazil[END_REF]. Moreover, the choice of the insulation thickness is a difficult task since it is a delicate compromise between the cost of insulation materials, which increases with the thickness, and the energy cost for heating/cooling the building, that decreases with the thickness [13,21].

The determination of the optimum thickness is a difficult task. It requires the computation of the heating or cooling loads of a building along a whole year by solving the heat transfer equation through building walls and roofs. Several studies reported in literature enhance the importance of efficient modeling to assess wall thermal loads. The influence of important wall parameters on the energy efficiency has been investigated: insulation properties in [START_REF] Ioannou | Energy performance and comfort in residential buildings: Sensitivity for building parameters and occupancy[END_REF], thermal inertia in [START_REF] Aste | The influence of the external walls thermal inertia on the energy performance of well insulated buildings[END_REF]39] or thermal conductivity in [START_REF] Kontoleon | The influence of concrete density and conductivity on walls' thermal inertia parameters under a variety of masonry and insulation placements[END_REF]. In [25], the effect of wall moisture content on the thermal response of the wall has been analysed. Some studies aimed at calibrating the results of the numerical models with observations in order to improve the reliability of the physical model, as for instance in [START_REF] Kim | Building energy model calibration with schedules derived from electricity use data[END_REF][START_REF] Mustafaraj | Model calibration for building energy efficiency simulation[END_REF][START_REF] Yang | A model calibration framework for simultaneous multi-level building energy simulation[END_REF][START_REF] Yuan | A simultaneous calibration and parameter ranking method for building energy models[END_REF].

These studies highlight the importance of efficient numerical models to predict heat transfer through walls. However, the elaboration of such models is a difficult issue since the numerical complexity arises from several points [START_REF] Capizzi | Thermal transients simulations of a building by a dynamic model based on thermal-electrical analogy: Evaluation and implementation issue[END_REF][START_REF] Mendes | Numerical methods for diffusion phenomena in building physics[END_REF]. First, the numerical model has to deal with different time scales since the computation is carried out over one year, as illustrated in [START_REF] Branco | Predicted versus observed heat consumption of a low energy multifamily complex in Switzerland based on long-term experimental data[END_REF], while the characteristic time scale of the heat diffusion scales with 15 min or even less. In addition, the boundary conditions are defined according to non-periodic climatic data [36]. Third, the wall configuration involves several layers inducing variations of the thermal properties with space. These difficulties imply an important numerical effort of the numerical model to compute the optimal thermal insulation.

In literature, several works circumvent the numerical difficulties to reduce the computational efforts. For instance, in [START_REF] Capizzi | Thermal transients simulations of a building by a dynamic model based on thermal-electrical analogy: Evaluation and implementation issue[END_REF], a whole building model is replaced by lumped approaches, denoted as RC model to reduced the computation efforts. This approach is used in many algorithms to solve the partial differential equation of heat transfer as mentioned in [START_REF] Fraisse | Development of a simplified and accurate building model based on electrical analogy[END_REF][START_REF] Naveros | Order selection of thermal models by frequency analysis of measurements for building energy efficiency estimation[END_REF][START_REF] Roels | On site characterisation of the overall heat loss coefficient: Comparison of different assessment methods by a blind validation exercise on a round robin test box[END_REF]. For the study of optimal thermal insulation, studies also attempt to avoid the numerical difficulties. For instance, in [START_REF] Ekici | A study on the optimum insulation thicknesses of various types of external walls with respect to different materials, fuels and climate zones in Turkey[END_REF], the simple degree-day method is used to compute the building energy consumption in four cities of Turkey. In [START_REF] Axaopoulos | Optimum insulation thickness for external walls on different orientations considering the speed and direction of the wind[END_REF]13,21] simulations are performed only on the hottest and coldest days/months of the year. In [START_REF] Ozel | Effect of wall orientation on the optimum insulation thickness by using a dynamic method[END_REF], the outside boundary condition is assumed as periodic by repeating a daily cycle of air temperature and solar radiation. In [3,[START_REF] Gasparin | LOCIE, UMR 5271[END_REF], a representative day for each month of the year is used to compute the yearly loads. In [START_REF] Daouas | Impact of external longwave radiation on optimum insulation thickness in Tunisian building roofs based on a dynamic analytical model[END_REF], an analytical solution is proposed for simulations carried out only for the hottest and coldest days of the year. For the climate of Paris, [START_REF] Berger | An innovative method for the design of high energy performance building envelopes[END_REF] proposed a parametric solution of the transient heat transfer using a reduced order model to determine the annual loads and consequently the optimum insulation thickness.

In this way, we propose an innovative method to simulate the transient heat transfer to determine the optimum insulation thickness, which will be tested in Brazilian buildings context. Instead simulating the hottest and coolest months we simulate the whole year to get the transmission loads. The method used for simulations is called Quasi-Uniform Nonlinear Transformation (QUNT) which is adaptive regarding the spatial grid [START_REF] Khakimzyanov | Numerical solution of conservation laws on moving grids[END_REF]. The grid points of the spatial mesh move following the monitor function which indicates their placement. The adapted grid is obtained according to the gradients of the problem indicated by solving a parabolic problem. As shown in [START_REF] Khakimzyanov | On supraconvergence phenomenon for second order centered finite differences on non-uniform grids[END_REF], this approach upgrades the usual finite-differences method, with a semi-implicit scheme, to the fourth order by using the non-uniform grid [START_REF] Blom | Algorithm 731; a moving-grid interface for systems of one-dimensional time-dependent partial differential equations[END_REF]. Just by changing the distribution of nodes one can gain two extra orders of accuracy, which allows simulating with efficiency using less discrete points. This numerical model enables to compute the temperature distribution for the entire year at a low computational effort. From these results, the optimal thermal insulation can be easily computed.

Therefore, this paper first describes in Section 2 the employed physical model of heat transfer. In Section 3, basics of the moving grid method is detailed before exploring the features of the scheme applied to a nonlinear case, presented in Section 4. Then, in Section 5, the method is applied to investigate the optimum insulation thickness in Brazilian buildings. The main conclusions of this study are outlined in Section 6.

Mathematical formulation of the physical problem

The physical problem considers one-dimensional heat transfer through a porous material. The spatial domain is defined by x ∈ [ 0, l ], in which x = 0 corresponds to the surface in contact with the inside room and, x = l, corresponds to the outside surface. The heat transfer is governed by conduction, which formulation is derived from Fourier's law:

c ∂T ∂t = ∂ ∂x k ∂T ∂x , (2.1)
where x is the spatial coordinate [m], t is the time [s], T is the temperature [K] * , k, the thermal conductivity [W/(m • K)] and c = ρ c p , the thermal storage, in which c p is the material specific heat [J/(kg • K)] and ρ, the material density [kg/m 3 ].

The Dirichlet-type boundary conditions are written as:

T (x = 0 , t) = T L (t) and T (x = l , t) = T R (t) .
Moreover, the internal Robin-type boundary condition is expressed as:

-k ∂T ∂x x = l = h R T -T R (t) ,
and on the external side, the boundary condition includes the convection, short-and longwave radiation:

k ∂T ∂x x = 0 = h L T -T L (t) -α q ∞ (t) + ǫ σ T 4 -T 4 sky (t) .
where T L and T R are the temperature of the air [K] that varies over time, h is the convective heat transfer coefficient [W/(m 2 • K)] and subscripts L and R represent the left and right boundaries. If the bounding surface is in contact with the outside building air, q ∞ is the total solar radiation [W/m 2 ] , which includes the direct, diffuse and reflexive radiations and α is the surface absorptivity [-]. For the long-wave radiation, T sky is the sky temperature

[K], ǫ is the emissivity [-] and σ is the Stefan-Boltzmann constant [W/(m 2 • K 4 )] .
For the initial condition, a linear temperature distribution in function of x is considered:

T (x , t = 0) = T i (x) . (2.
2)

The heat flux density q [W/m 2 ] at the point x j is computed as:

q j (t) def := -k ∂T ∂x (x , t) x = x j . (2.3)
Finally, the integration of the resulting heat flux density during a time period will give the transmission loads E [J/m 2 ] [8, 9]:

E = ˆt1 t 0 q (τ ) dτ , (2.4)
which can be evaluated on daily or monthly periods.

Method of solution

The problem previously described is solved by using the Quasi-Uniform Nonlinear Transformation (QUNT) [START_REF] Khakimzyanov | Numerical solution of conservation laws on moving grids[END_REF]. However, the problem must be expressed in a dimensionless formulation before being solved numerically. For example, this process can reduce the number of parameters of the problem allowing direct comparisons of models. For a detailed explanation of why this practice is important one can consult [START_REF] Gasparin | Stable explicit schemes for simulation of nonlinear moisture transfer in porous materials[END_REF].

In this way, we define the following dimensionless quantities:

u def := T T 0 , x ⋆ def := x l , t ⋆ def := t t 0 , c ⋆ def := c c 0 , k ⋆ def := k k 0 , Bi def := h • l k 0 , Fo def := t 0 • k 0 l 2 • c 0 , q ⋆ def := l • q T 0 • k 0 , R ⋆ lw def := ǫ • σ • l • T 3 0 k 0 ,
where the subscript 0 represents a reference value, chosen according to the application problem and the superscript ⋆ represents a dimensionless quantity of the same variable. Therefore, the governing Equation (2.1) can be written in a dimensionless form as:

c ⋆ ∂u ∂t ⋆ = Fo ∂ ∂x ⋆ k ⋆ ∂u ∂x ⋆ , (3.1)
with the following dimensionless formulation of the Dirichlet and Robin-type boundary conditions:

u (x ⋆ = 0 , t ⋆ ) = u L (t ⋆ ) , u (x ⋆ = 1 , t ⋆ ) = u R (t ⋆ ) , k ⋆ ∂u ∂x ⋆ x ⋆ = 0 = Bi L u -u L (t ⋆ ) -α q ⋆ ∞ (t ⋆ ) + R ⋆ lw u 4 -u 4 sky (t ⋆ ) , -k ⋆ ∂u ∂x ⋆ x ⋆ = 1 = Bi R u -u R (t ⋆ ) .
For the initial condition the dimensionless form is:

u (x ⋆ , t ⋆ = 0) = u i (x ⋆ ) .
Moreover, the dimensionless form of the heat flux density is written as:

q ⋆ j (t ⋆ ) = -k ⋆ ∂u ∂x ⋆ x ⋆ = x ⋆ j .
In the following, we drop ⋆ for the sake of notation compactness but without losing generality.

Quasi-Uniform Nonlinear Transformation

The Quasi-Uniform Nonlinear Transformation (QUNT) [START_REF] Blom | Algorithm 731; a moving-grid interface for systems of one-dimensional time-dependent partial differential equations[END_REF][START_REF] Khakimzyanov | On supraconvergence phenomenon for second order centered finite differences on non-uniform grids[END_REF][START_REF] Khakimzyanov | Numerical solution of conservation laws on moving grids[END_REF] is a numerical method with an adaptive spatial grid which provides a high resolution where it is required. The adaptivity is obtained by using a monitor function which indicate the location of the spatial grid points during the simulation. The adapted grid is obtained by means of the equidistribution principle. This moving grid nodes approach has the advantage of being conservative in space and also second-order accurate.

Before solving the problem on a moving mesh, one have to chose a robust and accurate scheme on a fixed grid. This scheme will be then generalized to incorporate the motion of the mesh. This scheme is described in the sequence.

A finite-difference scheme on a fixed uniform mesh

Consider the dimensionless form of the heat diffusion Equation (3.1) with a simplified notation, just for describing the numerical method:

c (u) ∂u ∂t = ∂ ∂x k (u) ∂u ∂x . (3.2)
The spatial domain I = [ 0 , l ] is uniformly discretized, with x j = j∆x, j ∈ Z being the spatial nodes, in which ∆x > 0 is the spatial step size. For the time domain P = [ 0, τ ], also consider a uniform discretization, where t n = n ∆t , n ∈ Z + are the temporal nodes and ∆t > 0 being the time step. Thus, by using an IMplicit-EXplicit (IMEX) scheme, Equation (3.2) is written in a discrete form:

c n j u n+1 j -u n j ∆t = 1 ∆x k n j+ 1 2 u n+1 j+1 -u n+1 j ∆x -k n j-1 2 u n+1 j -u n+1 j-1 ∆x .
This scheme approximates the continuous operator to order O (∆x 2 + ∆t). The advantage of this semi-implicit scheme over the fully implicit is to avoid sub-iterations in the solution procedure and, at the same time, being stable and consistent. For simplicity, the boundary conditions are not treated here. In what follows, the method based on this scheme with the moving grid is describe.

Finite-difference scheme on a moving grid

Consider the computational domain I = [ 0 , l ] , the reference domain Q = [ 0 , 1 ] and a bijective time-depend mapping from Q to I :

x (q, t) : Q × R + -→ I .
This represents the quasi-uniform nonlinear transformation. In addition, the boundary points are required to map into each other:

x (0 , t) = 0 and x (1 , t) = l .
The reference domain Q is uniformly discretized into N elements, with q j = j h, j ∈ Z + representing the nodes of the grid, which is equally spaced h = 1/N . However, only the image of nodes q j under the map x (q , t) are needed, since they constitute the nodes of the moving mesh:

x (q j , t n ) = x n j .
The diffusion Equation (3.2) is rewritten on the domain Q, with the help of the composed function v (q , t) def := (u • x) (q , t) ≡ u x (q , t) , t :

J (q , t) c (v) ∂v ∂t = ∂ ∂q k(v) J (q , t) ∂v ∂q + c (v) ∂v ∂q ∂x ∂t ,
where

J def := ∂x ∂q
is the Jacobian of the transformation x (q , t) .

Therefore, the discrete form of the diffusion heat equation on a moving mesh grid is:

J n j c n j v n+1 j -v n j ∆t = 1 h k n j+ 1 2 J n j+ 1 2 v n+1 j+1 -v n+1 j h - k n j-1 2 J n j-1 2 v n+1 j -v n+1 j-1 h + + c n j v n+1 j+1 -v n+1 j-1 2 h • x n+1 j -x n j ∆t .
In summary, this is the parametrization process of the mesh motion by a bijective mapping.

In the sequence, the construction of this map is explained.

Construction and motion of the grid

The equidistribution method is used to construct the adaptive grid. To control the distribution of the nodes, a positive valued function w (x, t) has to be chosen. This function is called monitor function. The choice of the monitor function and its parameters are important for the accuracy of the scheme. In particular, for unsteady nonlinear problems, there are several choice options. One can find more information, for example in [START_REF] Stockie | A Moving Mesh Method for Onedimensional Hyperbolic Conservation Laws[END_REF]. In this work, the function chosen is based on [START_REF] Khakimzyanov | Numerical solution of conservation laws on moving grids[END_REF] and can be written as:

w (x, t) = 1 + α 1 | u | β 1 + α 2 | ∂u ∂x | β 2 , (3.3)
where α 1 , α 2 , β 1 and β 2 are positive real parameters. The quantity w (x , t) assumes large positive values in areas where the solution has important gradients, inducing the elements of the grid mesh to approach. The discrete formulation of the monitor function presented in Eq. (3.3) is written as:

w n j = 1 + α 1 | u n j | β 1 + α 2 | u n j+1 -u n j x n j+1 -x n j | β 2 .
A non-uniform grid I is given if we construct the mapping x (q , t) : Q -→ I and evaluate it in the nodes of the uniform grid. To this end, the equidistribution method proposes that the desired mapping x (q , t) be obtained as a solution of a nonlinear parabolic problem, which is written as:

∂ ∂q w (x , t) ∂x ∂q = β ∂x ∂t , (3.4)
with x (0, t) = 0 and x (1, t) = l as the boundary conditions. Eq. (3.4) can be written in the discrete form as:

1 h w n j+ 1 2 x n+1 j+1 -x n+1 j h -w n j-1 2 x n+1 j+1 -x n+1 j h = β x n+1 j -x n j ∆t ,
with the same boundary conditions x n + 1 0 = 0 and x n + 1 N = l . The parameter β > 0 plays the role of the inverse diffusion coefficient and it controls the smoothness of nodes trajectory.

Initial grid generation.

As the problem depends also on time, it is of capital importance to obtain a high-quality initial mesh. The initial condition must be adapted to the new grid before starting the dynamical simulation. At t = 0 we compute the monitor function w (x, 0) of the initial condition u (x, 0) . Then, the mapping x (q, 0) is determined as the solution of Eq. (3.4), which is reduced to a simple second-order ordinary differential equation:

∂ ∂q w (x , 0) ∂x ∂q = 0 ,
supplemented with the Dirichlet boundary condition. Thus, the finite-difference approximation of this latter differential equation is:

1 h w 0 j+ 1 2 x 0 j+1 -x 0 j h -w 0 j-1 2 x 0 j+1 -x 0 j h = 0 ,
with the discrete boundary conditions x 0 0 = 0 , x 0 N = l . This latter nonlinear system of equations is solved iteratively, with iterations initialized with a uniform grid as the first guess. Its solution satisfies the equidistribution principle: in areas where w 0 j + 1 2 takes large values, the space between two neighboring nodes x 0 j+1 and x 0 j has to be inversely proportionally small. One can find more information in [START_REF] Khakimzyanov | On supraconvergence phenomenon for second order centered finite differences on non-uniform grids[END_REF].

Smoothing step. To ensure the smoothness of the mesh motion the monitor function w can be filtered. This process, called smoothing step enable to enlarge the set of acceptable values of the parameters α and β of the monitor function. According to [START_REF] Khakimzyanov | Numerical solution of conservation laws on moving grids[END_REF], the following implicit scheme can produce robust results:

w j + 1 2 = w j + 1 2 -σ w j + 1 2 + σ 2 w j -1 2 + w j + 3 2 , (3.5)
where σ is a positive smoothing parameter. To complete Eq. (3.5), boundary conditions are taken as in the original problem:

w 1 2 = w 1 2 , w N -1 2 = w N -1 2 .
The smoothed discrete monitor function

w j + 1 2 N -1 j = 0
is then obtained by solving the linear System (3.5). Note that the smoothing operator is applied to the monitor function.

On the choice of the parameters. The parameters α 1 , α 2 , β 1 , β 2 , β and σ are chosen according to the problem and the solution under consideration. To determine the optimal values one has to make some numerical experiments. Interested readers are invited to consult [START_REF] Khakimzyanov | Numerical solution of conservation laws on moving grids[END_REF] for more details on this aspect.

Validation of the numerical solution

All the numerical results in this paper were computed using Matlab™ [28]. The reference solution is computed using the open source package Chebfun [START_REF] Driscoll | Chebfun Guide[END_REF]. Using pde15s function, it enables to compute a numerical solution of a partial derivative equation with the Chebyshev polynomials adaptive spectral methods.

The Euclidean error (distance) for a determined field u at time t n is computed as:

δ (x) def := u ref (x , t n ) -u num (x , t n ) 2 ,
where u ref is the reference solution computed with Chebfun and u num is the solution computed with the QUNT or classical schemes. The error over the time is computed as a function of x by the following formulation:

ε (x) def := 1 N t N t n = 1 u ref (x , t n ) -u num (x , t n ) 2 ,
where N t is the number of temporal steps. The global L ∞ error is given by the maximum value of ε (x) :

ε ∞ def := sup x ∈ [ 0 , l ] ε (x) .
For the heat flux density, the Euclidean error (distance) is computed as:

ξ (t) def := q ref (t) -q num (t) 2 ,
where q ref is the reference flux density computed with Chebfun solution and q num is the flux density computed with the QUNT or classical solutions. The global L ∞ error of the flux density is also given by the maximum value but of the ξ (t) :

ξ ∞ def := sup t ∈ [ 0 , τ ] ξ (t) .

Numerical Benchmarking

The purpose of this section is to validate the QUNT method by comparing it with the IMEX scheme and the reference solution. The problem is inspired by [START_REF] Gasparin | Stable explicit schemes for simulation of nonlinear moisture transfer in porous materials[END_REF]. As the purpose is to evaluate the efficiency of the QUNT approach, and not analyzing the physical phenomena, the study is performed using the dimensionless problem described by Equation (3.1) with Dirichilet-type boundary conditions and the following dimensionless material properties:

k ⋆ (u) = 1 + 0.91 u + 600 • exp -10 u -1.5 2 , c ⋆ (u) = 900 -656 u + 10 4 • exp -5 u -1.5 2 .
The dimensionless ambient temperature at the left and right boundaries are:

u L (t ⋆ ) = 1 -0.5 • sin 2 π t ⋆ 12 and u R (t ⋆ ) = 1 + 0.5 • sin 2 π t ⋆ 24 ,
which are graphically given in Figure 1. The initial temperature is considered uniform over the material and its dimensionless formulation is given by:

u ( x, 0 ) = 1 .
The monitoring function is defined by Equation (3.3) with following values for its parameters:

α 1 = 0.9 , β 1 = 2 , α 2 = 0.1 , β 2 = 2 .
These parameters were determined by our numerical investigations and correspond to the best combination for this case study. For the grid diffusion parameter, the best value tested was β = 100 and for the smoothing parameter σ = 10 . Simulations are performed up to the final time of τ = 48 , with ∆t ⋆ = 5 • 10 -3 and a number of spatial points N x = 51 , at first fixed for both methods.

Results and discussion

The sample trajectory of the spatial nodes can be observed in Figure 2. The spatial points concentrate where the values of the gradient are high, particularly on the boundaries where important variations of the fields are imposed.

Figure 3(a) presents the temperature profile for t ⋆ = 16 of the computed solutions. As the material has a low thermal conductivity, the temperature in the interior of the material does not change as fast as it does at near the boundaries. Consequently, high gradients can be observed on the boundaries, making the spatial nodes concentrate in these areas. The QUNT method handles much better these high gradients than the IMEX, providing more accurate results. As the QUNT method was built based on the IMEX scheme, we can conclude that the moving grid modifications improved the numerical scheme. The absolute error δ for this instant of time is presented in Figure 3(b). This figure shows that the QUNT has a spatial accuracy the order of O (10 -4 ) , while the IMEX scheme has an accuracy the order of O (10 -3 ) . Thus for the same quantity of spatial nodes (N x = 51) , the QUNT can be 10 times more accurate than the numerical solution of the IMEX scheme. If we look at the boundaries, this ratio is even higher reaching ≈ 100 times. The results for the error over the time are presented in Figure 4. The order of the error remains stable for the whole simulation, with a maximum error of ε ∞ = 4.78 • 10 -4 for the QUNT and ε ∞ = 2.82 • 10 -2 for the IMEX method, which confirms ≈ 100 ratio between two approaches.

After simulating the temperature field, it is of great importance to get the values of the heat flux density at the boundaries, to compute, for instance, the transmission loads. Therefore, in what follows, we study the influence of the number of spatial nodes on the heat flux density.

Figures 5(a) and 5(b) present the heat flux evaluated at the left and right boundaries, respectively, with the derivative approximated with a finite-difference scheme second-order accurate in space. The results show that there is a difference between the flux computed for the uniform grid with the one computed for the non-uniform grid. Figures 6(a) and 6(b) represent the error of the fluxes ξ , computed with the QUNT and IMEX solutions for the left and right boundaries, respectively. In general, the flux computed with the non-uniform grid is more accurate than the flux computed with the uniform grid. By using 51 spatial nodes, the maximum error of the flux computed with an uniform grid is ξ ∞ = 1.9 • 10 -1 and the flux computed with a non-uniform grid is ξ ∞ = 8.5 • 10 -2 .

Convergence study

A convergence study was carried out to compare the QUNT with the IMEX scheme and to analyze the accuracy of both temperature distribution field and heat flux density. Thus, for different values of ∆t ⋆ = 10 -1 ; 10 -2 ; 10 -3 ; 10 -4 , the global error ε ∞ was computed as a function of the number of spatial nodes N x ∈ [ 10 , 80 ]. It can be seen in Figure 7 that for N x < 20 both methods provide solutions of the field with approximately the same order of accuracy, which is insufficient to get good approximations of the solution. For the IMEX method, less than 100 nodes are not enough to give a solution with an accuracy better than O (10 -2 ) . Even if the ∆t ⋆ is decreased, the global error is still of the same order of accuracy due to insufficient spatial discretization. On the other hand, for the QUNT method, as the number of spatial nodes increase, the global error decreases faster than the IMEX method. For the QUNT adaptive grid, with only 40 nodes, an accuracy of O (10 -3 ) can be achieved for the three values of ∆t ⋆ .

Figures 8(a), 8(b), 8(c) and 8(d) display the results for the global error of the flux density ξ ∞ , as a function of the number of spatial nodes N x ∈ [ 20 , 100 ] also for different values of ∆t ⋆ = { 10 -1 ; 10 -2 ; 10 -3 ; 10 -4 }. Results show that the computation of the flux density is less accurate than the computation of the temperature field. This happens since the error increases when the derivative has to be computed. The QUNT method also provides solutions for the flux more accurate than the IMEX method. As the analysis of the error for the field, the QUNT approach proved to be more sensible regarding variations of ∆t, as the solution becomes smoother and the moving grid adapts better when ∆t decreases. Although, the same accuracy is observed for ∆t = 10 -3 and ∆t = 10 -4 because the order of the error is also related with the spatial discretization. To obtain higher accuracy, it is necessary to increase the derivative approximation order. In this case, the error of the flux is relatively high due to the nonlinearity of the case study. For a less nonlinear case, the error of the flux density can be considerably reduced.

Therefore, it has been clearly noticed the QUNT approach computes an accurate solution of the problem. Moreover, the moving grid method converges to an accurate solution with less spatial nodes than the method with an uniform grid. In a long-term simulation, it does provide significant gains in regard to the CPU usage. The computer used to assess the computational cost has a processor Intel ® Core TM i5 @ 2.80GHz×4 with 15.6GB of RAM.

To evaluate the gains in terms of computational cost, Table 1 provides the CPU time to compute the solution using the Crank-Nicolson, the QUNT and the IMEX methods, which have been evaluated using the Matlab™ platform. All methods computed a solution with the same order of accuracy ε ∞ ∼ = O (10 -4 ) . For this, it was used a time step of ∆t ⋆ = 10 -2 , and 51 spatial nodes for the QUNT and 1 001 and 501 spatial nodes for the Crank-Nicolson and IMEX methods, respectively. The computational effort to perform the simulation increases linearly with the simulation time. The QUNT and IMEX methods can compute the solution 25% faster than the Crank-Nicolson approach. It is preferable to focus on the ratio of computer run-time rather than on absolute values, since it is system-dependent. The Crank-Nicolson scheme was implemented with a fixed point approach to treat the non-linearity and it was used here only for comparison purposes. One should notice that the codes are not optimized, and the differences between the CPU time of QUNT and IMEX methods are low, even with the QUNT method performing two extra matrix inversions.

Optimum insulation in Brazil

In the previous section, the QUNT method was presented and validated. Now, the method is used to compute the solution of a transient one-dimensional case study of heat transfer through building walls for a period of one year. The results are obtained with less computational efforts and enable to estimate the optimum thermal insulation for large cities in Brazil.

City

T min Climatic zones. Brazil has a total surface area of 8.5 million km 2 , the fifth largest country in the world, occupying almost half of the entire South America continent. In such a large country a variety of climates exist. However, most of the country can be defined as being tropical and sub-tropical and energy consumption for space conditioning is dominated by cooling. Although, in the Southern region of Brazil, energy is also used for heating during the winter months (June -August).

[ • C] T mean [ • C] T max [ • C] Altitude [m] Location Clim. zone Curitiba - 2
For this study, four Brazilian cities -Curitiba, Rio de Janeiro, São Paulo and Salvador -have been selected for determining their optimal insulation thickness, chosen due to the diversity of the climate and due to their large population [START_REF] Morishita | Issues about moisture in residential buildings of Brazil[END_REF]. However, the method used in this study can be applied to any city desired. Based on ASHRAE Standard 90.1 -2013 [START_REF]Chapter 26: Heat, Air, and Moisture Control in Building Assemblies -Material Properties[END_REF], these cities are located in three different climatic zones: Very Hot-Humid (zone 1), Hot-Humid (zone 2) and Warm-Humid (zone 3). Table 2 presents some characteristics of the chosen cities.

Simulation parameters. Bricks are commonly used in typical wall structures in Brazilian constructions [START_REF]Desempenho térmico de edificações Parte 1: Definições, símbolos e unidades[END_REF][START_REF] Morishita | Issues about moisture in residential buildings of Brazil[END_REF]. As illustrated in Figure 9, we consider 3 types of walls: Wall-1 is composed only of brick; Wall-2 includes also a layer of insulation on the internal side; Wall-3 has a layer of insulation on the external side. The brick is 15 cm thick, while the thickness of insulation material (extruded polystyrene) varies. The thermal properties of the materials are gathered from the IEA Annex- 24 [27] and are reported in Table 3. As the initial condition, a steady-state temperature distribution within the wall is assumed:

T i (x) = T R (0) -T L (0) l x + T L (0) [ • C] .
Climatic conditions of temperature and solar radiation are considered on the outside part of the wall for a period of one year. The outside temperature and total radiation are given in Figures 10 and11. The outside heat flux is computed using radiation data from Domus [START_REF] Mendes | A simulation environment for performance analysis of HVAC systems[END_REF] for a South-oriented wall. The solar absorptivity of the outside surface is taken as α = 0.5 . The indoor air temperature T R (t) has a sinusoidal variation as illustrated in Figure 10. During the summer period, the temperature is set to 25 • C and during the winter to 20 • C [START_REF] Melo | Building energy performance assessment: Comparison between ASHRAE standard 90.1 and Brazilian regulation[END_REF]. The convective heat transfer coefficients at the outer and the inner surfaces are set to constant values h L = 25 W/(m 2 • K) and h R = 10 W/(m 2 • K), respectively. The long-wave radiative transfer on the external surface is only considered for the roof simulation.

According to Figure 10, the inside temperature in all cities varies in the same way. What can be observed, is that Curitiba is the coldest city, which means it needs more heating than cooling. In Rio de Janeiro and Salvador, the external temperatures are higher than the inside temperature, which means that the cooling system predominates. In addition, Salvador has a higher mean external temperature than Rio de Janeiro. São Paulo has the mildest weather among the four cities, but the external temperatures are mostly lower than the internal ones, which makes the heating need prevails.

In next sections, the thermal behavior of buildings located in these cities is studied for different wall configurations, various thicknesses of insulation material and different wall orientations, with respect to the heat flux density q and transmission loads E on the inside surface of the building wall. Three cases are distinguished. First, the benefits of insulating buildings are demonstrated. The second case analyzes if the insulation should be placed on the inside or outside part of the building. Last, the effect of the wall orientation is studied.

First Case: insulated vs. non-insulated buildings

For the first simulation, two South-oriented walls in the city of Curitiba are considered. One has no insulation (Wall-1) and the other one has 10 cm insulation (Wall-2). Climatic conditions of the temperature and solar radiation are given in Figures 10(a With the numerical values provided previously, Problem (2.1) is solved using the QUNT method. The monitoring function is defined by Equation (3.3), which has the following values for its parameters: α 1 = 0.8 , β 1 = 2 , α 2 = 0.2 , β 2 = 3 . The grid diffusion and the smoothing parameters are, respectively, β = 50 and σ = 5 . The solution is obtained for a number of N x = 41 spatial nodes with a time discretization of ∆t ⋆ = 0.1 , which is equivalent to ∆t = 6 min.

The temperature is computed for all points in the single-and multi-layered wall configurations. Then, the heat flux on the inside surface of the building wall is estimated and presented in Figure 12(a). For the case with no insulation, the heat flux has an important variation, alternating between positive and negative values. With 10 cm insulation, it strongly modifies the thermal behavior of the wall and the heat flux is considerably reduced in this case.

Yearly transmission loads are calculated based on Equation (2.4). Hourly variation of inside surface heat flux density is integrated over a 24 h period to obtain a daily total load. Then, daily total loads are integrated to get the monthly heating/cooling transmission loads presented in Figure 12(b), for Wall-1 and Wall-2 . The heating/cooling loads are ten times higher for the wall without insulation indicating that much more energy is spent in order to maintain the temperature inside of the building.

With these results, the first conclusion is that insulation may play an important role in the construction sector in Brazil. In what follows, two configurations are tested, one with the insulation inside of the building wall (Wall-2) and another one with the insulation outside (Wall-3).

Second Case: inside insulation vs. outside insulation

For this second case, it is also considered a South-oriented wall in the city of Curitiba, with 10 cm insulation. Two wall configurations are verified: one with the insulation inside the building (Wall-2) and another one with the insulation outside (Wall-3). The same climatic conditions and numerical values of the previous case (Section 5.1) are taken.

Figure 13(a) shows the effect of the insulation position on the heat flux density during one year of simulation. The wall with external insulation has the lowest heat flux values but with a small discrepancy, with the mean difference between the heat flux of Wall-2 and Wall-3 being of the order of 0.1 W/m 2 , for a 10 cm insulation. This difference can also be noted in the monthly heating/cooling transmission loads as shown in Figure 13(b). Depending on the month, the value of the energy load can be higher for the insulation outside of the building wall. However, the annual heating/cooling transmission loads are very close, with 33.83 MJ/m 2 for Wall-2 and 34.81 MJ/m 2 for Wall-3 . In [START_REF] Axaopoulos | Optimum insulation thickness for external walls on different orientations considering the speed and direction of the wind[END_REF] they also have similar results for the annual loads regarding the position of the insulation (inside/outside), even with different assumptions and climatic conditions.

J F M A M J J A S O N D 12 -5 -4 -3 -2 -1 0 1 2 3 (a)

Parametric study

To improve this analysis, a parametric study is performed, in which the thickness of the insulation material varies on the interval l i = [0.01, 0.3] m , for both wall configurations (Wall-2 and Wall-3).

Simulations are performed with the same numerical and physical configuration of the previous case studies. The annual heating/cooling transmission loads are calculated and presented as a function of l i in Figure 14. As expected the thicker the insulation, the lower the energy consumption for either heating or cooling. One can easily notice that transmission loads drop fast for small insulation thickness values. Besides, for thin insulations, the difference between the two configurations is more evident. Although, this difference is minimized for high values of thicknesses.

With those results, another conclusion from this case study is that the difference of configurations with the insulation inside or outside does not cause important changes on annual heating/cooling when the thickness is higher than 2 cm, despite the thermal mass on the internal side provides important effects on a shorter timescale. For this reason, just for simulation purposes, in what follows, the configuration of Wall-2 is used to study the impact of wall-orientation.

Third Case: wall orientation influence

For the third case study, the same features of the previous case are considered but changing only the facing wall orientation (North, South, East, and West). As Brazil is located in the Southern hemisphere, the North façade wall of a building is the one that receives more short-wave radiation.

Figure 15 shows the effect of wall orientation on the monthly heating/cooling transmission loads of Wall-2 model with a 10 cm thick insulation. The monthly transmission loads values are mostly negative indicating a heating demand, which is compatible with the weather of Curitiba. In colder months, the South façade has the highest heating transmission loads because of the low solar heat gain, which is the opposite for the North-facing wall.

The effect of wall orientation on yearly heating and cooling transmission loads are presented in function of the insulation thickness in Figure 16. The South-facing wall has the highest transmission loads, then, East-, North-and West-facing walls with very similar values. As previously observed, transmission loads quickly drop is faster for small values of insulation thickness.

J F M A M J J A S O N D

The roof

The roof is one of the building components whereby occurs the most significant heat gain/loss. An identical analysis of the wall structure is performed for a roofing system. For this analysis, a flat roof is taken into account and we consider configurations identical to the wall structure, with internal and external insulation. The roofing is made of a 15 cm concrete fixed layer while the thickness of the insulation material varies. The thermal properties of the materials are reported in Table 3.

Simulations were performed by considering radiative and convective heat transfer at the external roof surface, which exchanges heat with the air and the sky, and receives solar radiation. The solar absorptivity of the outside surface of the roof is taken as 0.5 . For the long-wave radiation exchange, a thermal emissivity of 0.9 is considered and the sky temperature is provided by the software Domus [START_REF] Mendes | A simulation environment for performance analysis of HVAC systems[END_REF]. Instead of computing the transmissions loads only for Curitiba, the four cities were comprised in this study.

Figure 17 presents the annual heating/cooling transmission loads as a function of insulation thickness, for the four Brazilian cities. Salvador and Rio de Janeiro are the cities that have the highest transmission loads, which means that they will need a thicker insulation. The differences for the transmission loads are higher for values of insulation thickness. It can also be noticed that the differences are reduced, while the insulation thickness increases.

Optimal thermal insulation: economic analysis

To assess the economic viability of the insulation, a simple economic analysis based on [START_REF] Berger | An innovative method for the design of high energy performance building envelopes[END_REF][START_REF] Wang | Handbook of Air Conditioning and Refrigeration[END_REF] is performed. The total cost of the insulation C T is estimated as:

C T = C E + C I [$/m 2 ] ,
where C E is the annual cost of the electric energy consumed and C I is the the cost of the insulation. The energy cost is calculated by using the annual conduction load E , the energy efficiency η and the energy price p e :

C E = E • p e η .
The insulation cost p i comprises the price of the material as well as its installation and it varies according to the thickness l i : The optimal insulation thickness is given when C T has minimum owning and operating costs. For this analysis, the parameters used are given in Table 4. The energy price is the average price of the electric power in Brazil for 2017 [5].

C I = p i • l i .
Figure 18 presents the insulation cost, the energy consumption cost, and the total cost, as a function of the insulation thickness. As the insulation thickness increases, the energy cost decreases because less energy is used to maintain the temperature inside of the building. However, installation cost increases linearly with the insulation thickness. Thus, the compromise between the energy cost and the insulation cost is the minimum value of their sum. For the case study of Section 5.2.1, it corresponds to an optimum insulation thickness of 0.5 m for both wall configurations (Wall-2 and Wall-3).

This economy study is extended to other three cities, considering the four oriented façades and with the insulation on the inside part of the building wall. the values of optimum insulation thickness for East, West, South, and North oriented walls for the four Brazilian cities. In the case of the wall, the thickness of the insulation is the same for the external and internal configurations. Results show that Salvador is the city that needs the largest insulation thickness, as the inside temperature has the largest difference regarding the outside temperature. Salvador is the hottest city among them with cooling demand along the whole year. The Northern and Western façades of Rio de Janeiro and São Paulo are the ones receiving the most solar radiation. In Rio de Janeiro, due to the cooling demand, these façades require larger insulation thickness, as reported in For the roof, the position of the insulation generate differences in the optimum thickness. As shown in Table 6, Curitiba and São Paulo have different values of the optimal thickness of external and internal roof insulation. For Rio de Janeiro and Salvador, the thickness of the insulation must be higher, as they exchange more heat through the roof. For all cities, the roof needs more insulation than the wall, as previously mentioned.

Conclusion

In this work, the Quasi-Uniform Nonlinear Transformation (QUNT) method has been presented as an innovative method to perform the transient heat transfer, which was applied to an energy optimization problem. The method is based on a non-uniform adaptive grid technique that identifies were the spatial nodes must be placed. It was first compared with a reference solution to verify its features. As a result, a satisfactory accuracy is shown by using much less spatial nodes than traditional methods, by moving its nodes where the temperature gradients are higher i.e., an average ratio of 100 compared with a uniform grid approach. Another remarkable point of this method is that it is easy to be implemented, by adding extra differential equations without being computational costly. An expressive gain of 25% is obtained.

The high accuracy of the solution instigated the analysis of the physical behavior. The optimum insulation thickness of buildings walls is determined by taking into account the wall orientation and the position of the insulation in Brazilian buildings, which is determined through a parametric study. The annual heating/cooling transmission loads have been calculated using weather data and transient heat flow, simulated over one whole year.

First, a case was focused on comparing a wall with insulation and another one with no insulation, applied to the climatic conditions of one Brazilian city. This case proved the importance of the insulation in reducing the transmission loads. Then, a second case study is simulated to verify if the insulation is better inside or outside of the building wall. The last case considered analyzing the effect of wall orientation and the roof assembly. As Brazil is mostly located in the Southern Hemisphere, the North façade is the one receiving more solar radiation. In consequence, for the hottest climate cities, the Northfacing wall is the one that requires the thicker insulation layers, while, in the South region -Curitiba -, the North-facing wall requires thinner layers of insulation.

A basic economic analysis was presented, demonstrating that the implementation of an insulation layer appears to be a cost-effective measure to save energy consumption. The optimum insulation thickness for all the wall and roof configurations studied is to be between 2 cm and 7 cm. The city that requires less insulation is São Paulo, due to its mild climate, and the city that needs more insulation is Salvador, due to its temperature to be high above the thermal comfort values.

As the proposed numerical method looked very promising, further work is recommended to integrate it to whole-building simulation tools so that more realistic cases, including internal gains and complex boundary conditions, can be taken into account in order to provide more conclusive results and enhance building energy simulation tools.

Nomenclature

Contents 1 1

 11 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Mathematical formulation of the physical problem . . . . . . . . . . . . . . . 3 Method of solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Quasi-Uniform Nonlinear Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A finite-difference scheme on a fixed uniform mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Finite-difference scheme on a moving grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Construction and motion of the grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Initial grid generation. 10 Smoothing step. 10 On the choice of the parameters. 11 3.2 Validation of the numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Numerical Benchmarking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Convergence study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Optimum insulation in Brazil . . . . . . . . . . . . . . . . . . . . . . . . . . . First Case: insulated vs. non-insulated buildings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Second Case: inside insulation vs. outside insulation . . . . . . . . . . . . . . . . . . . . . . . . . . Parametric study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Third Case: wall orientation influence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The roof . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Optimal thermal insulation: economic analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 .

 1 Figure 1. Dimensionless ambient temperatures.

Figure 2 .

 2 Figure 2. Trajectory of the grid nodes in space-time.

Figure 3 .

 3 Figure 3. Temperature profiles for the time t ⋆ = 16 (a) and the error δ (x ⋆ ) for the same profile (b).

Figure 4 .

 4 Figure 4. The error ε (x) over the whole simulation.

Figure 5 .Figure 6 .

 56 Figure 5. Flux evolution computed at the left (a) and at the right (b) boundaries.

  Figures 7(a), 7(b), 7(c), and 7(d) display the results.

Figure 7 .

 7 Figure 7. Global error of the temperature field ε ∞ computed in function of the number of spatial nodes, for ∆t ⋆ = 10 -1 (a), for ∆t ⋆ = 10 -2 (b), for ∆t ⋆ = 10 -3 (c), and for ∆t ⋆ = 10 -4 (d).

Figure 8 .

 8 Figure 8. Global error of the heat flux density ξ ∞ computed as a function of the number of spatial nodes, for ∆t ⋆ = 10 -1 (a), for ∆t ⋆ = 10 -2 (b), for ∆t ⋆ = 10 -3 (c) and for ∆t ⋆ = 10 -4 (d).

Figure 9 .

 9 Figure 9. Schematic representation of the wall's composition.

Material ρ 0

 0 [m 3 /kg] c p [J/kg/K] k [W/m/K]

  ) and 11.

Figure 10 .

 10 Figure 10. Inside and outside air temperatures for Curitiba (a), Rio de Janeiro (b), Salvador (c) and São Paulo (d).

Figure 11 .

 11 Figure 11. The heat flux incidence on the South wall façade.

Figure 12 .

 12 Figure 12. Heat flux density on the inside wall surface (a) and evolution of the monthly conduction load, for Wall-1 and Wall-2 models.

Figure 13 .

 13 Figure 13. Heat flux density on the inside wall surface (a) and monthly conduction loads (b) for Wall-2 and Wall-3.

Figure 14 .

 14 Figure 14. Annual heating and cooling transmission loads for Wall-2 and Wall-3 in function of insulation thicknesses.

Figure 15 .

 15 Figure 15. Monthly heating and cooling transmission loads for Curitiba at different orientations.

Figure 16 .

 16 Figure 16. Effect of wall orientation on annual cooling and heating loads in Curitiba for different insulation thickness.

Figure 17 .

 17 Figure 17. Annual heating/cooling loads associated to the roof for different insulation thicknesses.

Table 1 .

 1 Computational time required for the numerical schemes to perform the nonlinear case for ε ∞ ∼ = O (10 -4 ) .

	Simulation time Crank-Nicolson QUNT IMEX
	2 days	1.43 s	1.12 s 1.08 s
	1 month	21.36 s	16.2 s 15.2 s
	1 year	259.32 s	205.1 s 185.9 s
	Spatial nodes	1001	51	501
	CPU time (%)	100%	75%	72%

Table 2 .

 2 Climate and geographic characteristic of the Brazilian cities.

Table 3 .

 3 Thermophysical properties of the materials [27].

Table 4 .

 4 Economic parameters.

	Insulation price Heating/cooling system efficiency Energy price
	p i	η	p e
	100 $/m 3	0.8	0.218 $/(kWh)

Table 5

 5

	summarizes

Insulation cost, energy cost and total cost on function of insulation thickness for Curitiba city in a South-facing wall.

Table 5 .

 5 Optimum insulation thickness according to wall orientation for inside (or outside) insulation.

	Façade orientation Curitiba Rio de Janeiro Salvador São Paulo
	South	5 cm	4 cm	5 cm	3 cm
	North	4 cm	5 cm	6 cm	2 cm
	West	4 cm	5 cm	6 cm	2 cm
	East	4 cm	4 cm	6 cm	3 cm

Table 5 .

 5 However, in São

Table 6 .

 6 Optimum insulation thickness according to roof insulation position.Paulo the heating need prevails, and both façades orientations have smaller thickness values since their transmission loads are lower. On the other hand, in Curitiba and Salvador, the heat gains are higher for the North, West and East façades, making the South-facing wall require thicker insulation in Curitiba and thinner in Salvador.

	Insulation position Curitiba Rio de Janeiro Salvador São Paulo
	Outside	3 cm	6 cm	7 cm	2 cm
	Inside	4 cm	6 cm	7 cm	3 cm

Acknowledgements

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES) -in the framework of the International Cooperation Program CAPES/COFECUB (Grant #774/13). The Authors also acknowledge the support from CNRS/INSIS and Cellule Énergie under the grant MN4BAT-2017. Professor Mendes thanks the Laboratory LAMA UMR 5127 for the warm hospitality during his visits in 2018, which were supported by the project MN4BAT under the AAP Recherche 2018 programme of the University Savoie Mont Blanc. Finally, the authors acknowledge the Junior Chair Research program "Building performance assessment, evaluation and enhancement" from the University Savoie Mont Blanc in collaboration with The French Atomic and Alternative Energy Center (CEA) and Scientific and Technical Center for Buildings (CSTB).