© А О Бекетаева

P Bruel

А Ж Найманова

Детально-сравнительный анализ взаимодействия сверхзвукового потока с поперечной газовой струей при больших параметрах нерасчетности

Проведено исследование взаимодействия пространственного сверхзвукового турбулентного потока газа с вдуваемой перпендикулярно со стенки звуковой струей, направленное на выявление и углубленное понимание механизмов образования вихрей за вдуваемой звуковой струей в сверхзвуковом натекающем потоке в зависимости от (отношение давления в струе к давлению в потоке). Решение трехмерных осредненных по Фавру уравнений Навье-Стокса, замкнутых k-ω моделью турбулентности, осуществляется с помощью алгоритма, построенного на основе ENO-схемы высокого порядка аппроксимации. Показано присутствие известных из ряда теоретических работ вихревых структур. Определены параметры нерасчетности, при которых появляется дополнительная пара вихрей, где одна пара возникает в результате взаимодействия замедленного потока струи за диском Маха и высокоскоростного восходящего внешнего потока струи, а вторая обусловлена взаимодействием потока над струей с бочкообразной структурой. В результате сравнительного анализа выявлены условия, при которых наблюдается четкая картина дополнительных роговых вихрей вблизи стенки в области за струей. Получена зависимость угла наклона головного скачка уплотнения от параметра нерасчетности. Установлено удовлетворительное согласие распределения давления на стенке перед струей в плоскости симметрии с экспериментальными данными.

Ключевые слова: численное моделирование, сверхзвуковое течение, совершенный газ, пограничный слой, уравнения Навье-Стокса, параметр нерасчетности, ударная волна. DOI: 10.21883/0000000000

Введение

Течения высокоскоростных струй в поперечном потоке обеспечивают эффективное смешивание топлива и окислителя, имеющего решающее значение для сверхзвукового горения. На данный момент известно большое количество экспериментальных [START_REF] Глаголев | Панов Ю[END_REF][2][3][4][5] и теоретических [6][START_REF] De-Chaun | CAI Ti-min[END_REF][START_REF] Khali | Mixing flow characteristics for a transverse sonic jet injecting into a supersonic crossflow[END_REF][START_REF] Viti | [END_REF][START_REF] Lu | Topology of supersonic jet interaction flowfields at high pressure rasio // Presented at the July 1-4[END_REF][START_REF] Бекетаева | // Прикладная механика и техническая физика[END_REF][12][START_REF] Бекетаева | // Прикладная механика и техническая физика[END_REF][START_REF] Волков | // Прикладная механика и техническая физика[END_REF] исследований, в которых достаточно хорошо изучен механизм образования ударно-волновой структуры струйного взаимодействия для умеренных параметров нерасчетности n (отношение давления в струе к давлению в потоке). Общую картину течения схематически можно представить следующим образом (рис. 1, a -пространственная картина, b -на линии симметрии): линия 1 -головной скачок уплотнения, возникающий вследствие торможения набегающего потока перед струей, линии 2, 3 -косой и замыкающий скачки уплотнения соответственно. Головной, косой и замыкающий скачки уплотнения, пересекаясь в одной точке, образуют λ-образную структуру. Буквами D и B обозначены диск Маха и бочкообразная структура в струе соответственно. Линиями 4 и 5 изображены хорошо известные два противоположно вращающихся вихря перед струей, появляющихся из-за большого градиента давления, создаваемого на стенке перед вдуваемой струей и, как следствие, отрыва пограничного слоя. На рис. 1, c схематически иллюстрируется вихревая структура за струей, которая, как известно, вносит существенный вклад в смешение топлива и окислителя. Здесь линией 6 представлен вихревой след, 8это общеизвестная противоположно-вращающаяся пара вихрей в самой струе, образованная прогибом струи и конвекцией поперечного потока, также здесь приведен подковообразный вихрь 7.

Однако в работе [START_REF] Viti | [END_REF] в результате численного моделирования течения с параметром нерасчетности n = 282 и в расчетах [START_REF] Бекетаева | // Прикладная механика и техническая физика[END_REF], проведенных с параметрами нерасчетности в диапазоне 10 ≤ n ≤ 50, выявлено наличие дополнительных пар вихрей за струей, которая схематически также иллюстрируется на рис. 1, c. В этой картине новыми являются вихревые структуры 9 и 10, где пара вихрей 9 возникает за счет взаимодействия струи, прохо-1514 А.О. Бекетаева, P. Bruel, А.Ж. Найманова дящей через диск Маха с восходящим внешним потоком, а вихри 10 зарождаются за счет взаимодействия потока с бочкообразной структурой в струе.

Следует также отметить, что исследования авторов [START_REF] Lu | Topology of supersonic jet interaction flowfields at high pressure rasio // Presented at the July 1-4[END_REF] для параметра нерасчетности 5 ≤ n ≤ 2000 подтверждают наличие известных вихрей и показывают дополнительную пару вихрей (роговых) за вдуваемой струей непосредственно возле стенки, где пограничный слой обтекает струю, возникающую как в дозвуковом, так и в сверхзвуковом течениях, однако в работах [START_REF] Viti | [END_REF][START_REF] Бекетаева | // Прикладная механика и техническая физика[END_REF] эти вихри не были выявлены. Представляется затруднительным проведение сравнительного анализа полученных данных с результатами работы [START_REF] Lu | Topology of supersonic jet interaction flowfields at high pressure rasio // Presented at the July 1-4[END_REF], так как в этой работе не приведено количественное описание возникновения и распространения вихрей. Все это свидетельствует о необходимости дальнейшего исследования вихревых структур в широком диапазоне параметра нерасчетности.

Целью исследования настоящей работы является численное изучение и детально-сравнительный анализ механизмов формирования вихревых и ударно-волновых структур в широком диапазоне параметра нерасчетности, а также получение обобщенных зависимостей характеристик потока в зависимости от этого параметра. Для этого численно моделируется задача вдува звуковых струй из круглых отверстий, расположенных симметрично на верхней и нижней стенках канала, перпендикулярно сверхзвуковому потоку. Для удобства вычисления рассматривается вдув струи только с нижней стенки. Схема течения и качественная картина вихревой структуры показаны на рис. 1.

Постановка задачи

Исходной является система трехмерных осредненных по Фавру уравнений Навье-Стокса для сжимаемого турбулентного газа, записанная в декартовой системе координат в консервативной форме

∂U ∂t + ∂(E -E ν ) ∂x + ∂(F -F ν ) ∂z + ∂(G -G ν ) ∂y = S, (1)
компоненты векторов U, E, F, G определяются выражениями

U =           ρ ρu ρv ρw E t ρk ρω           , E =           ρu ρu 2 + P ρuv ρuw (E t + P)u ρuk ρuω           , F =           ρw ρuw ρvw ρw 2 + P (E t + P)w ρwk ρwω           , G =           ρv ρuv ρv 2 + P ρvw (E t + P)v ρvk ρvω           , а компоненты E v , F v , G v связаны с вязкими напряжени- ями E v = 0, τ x x , τ x y , τ x z , uτ x x + vτ x z + wτ x z -q x , 1 Re (µ l + σ k µ t ) ∂k ∂x , 1 Re (µ l + σ ω µ) ∂ω ∂x T , F v = 0, τ x z , τ y z , τ z z , uτ x z + vτ y z + wτ z z -q z , 1 Re (µ l + σ k µ t ) ∂k ∂z , 1 Re (µ l + σ ω µ) ∂ω ∂z T , G v = 0,
q x = µ (γ -1)M 2 ∞ Pr Re ∂T ∂x , q y = µ (γ -1)M 2 ∞ Pr Re ∂T ∂y , q z = µ (γ -1)M 2 ∞ Pr Re ∂T ∂z .
Здесь µ = µ l + µ t , где µ l -молекулярная вязкость (определяется по формуле Сазерленда), µ t -турбулентная вязкость. Вектор дополнительных членов имеет следующий вид: Здесь k, ω -кинетическая энергия турбулентности и скорость диссипации кинетической энергии турбулентности, P k -член генерации турбулентности, турбулентная вязкость определяется по формуле µ t = ρk ω [START_REF] Wilcox | A two-equation turbulence model for wallbounded and free-share flow[END_REF]. Выражения для давления и температуры запишутся следующим образом: Граничные и начальные условия имеют следующий вид: на входе, а также в качестве начальных данных задаются параметры потока

S = 0, 0, 0, 0, 0, (P k -β * ρωk), (γ * ρP k /µ t -βρω 2 ) T , где P k = µ t ∂u i ∂x j + ∂u j ∂x i ∂u i ∂x j - 2 3
∂u k ∂x k 2 - 2 3 ρk ∂u k ∂x k , i, j, k = 1, 2,
P = (γ -1) E t - 1 2 (ρu 2 + ρw 2 + ρv 2 ) , T = 1 ρc v E t - 1 2 (ρu 2 + ρw 2 + ρv 2 ) , c v = 1 γ(γ -1)M 2 ∞ . Здесь t -время, u, w, v -компоненты скорости потока в продольном и поперечных направлениях, ρ - плотность, E t -полная энергия, P -давление, T - температура, c v -теплоемкость при постоянном объ- еме, γ -показатель адиабаты, M ∞ -число Маха потока, далее индекс отнесен к параметрам струи, ∞ - к параметрам потока.
u = 1, v = 0, w = 0, ρ = 1, T ∞ = 1, x = 0, 0 ≤ y ≤ H y , 0 ≤ z ≤ H z ,
Начальные данные для k, ω параметров определялись с использованием алгебраической модели турбулентности Болдуина-Ломакса по известным осредненным физическим параметрам входного потока и подробно приведены в работе [START_REF] Бекетаева | // Прикладная механика и техническая физика[END_REF].

Вблизи стенки задается пограничный слой, толщина которого определяется по формуле δ 1 = 0.37x(Re) -0.2 . Также задается пристенный слой (10% от пограничного слоя) δ 2 = 0.1δ 1 . Продольная составляющая скорости u, а также температура и плотности аппроксимируются в соответствии с работой [START_REF] Бекетаева | // Прикладная механика и техническая физика[END_REF]:

А.О. Бекетаева, P. Bruel, А.Ж. Найманова на струе: u = 0, v = 0, w = T 0 M 0 /M ∞ , P 0 = nP ∞ , T 0 = 3.71, z = 0, |x 2 + y 2 | ≤ R. Здесь n = P 0 /P ∞ -параметр нерасчетности.
на верхней границе задается условие симметрии:

w = 0; ∂u ∂z = 0; ∂v ∂z = 0; ∂T ∂z = 0; ∂k ∂z = 0; ∂ω ∂z = 0; z = H z , 0 ≤ x ≤ H x , 0 ≤ y ≤ H y ;
на боковых границах:

∂u ∂y = ∂v ∂y = ∂w ∂y = ∂ρ ∂y = ∂k ∂y = ∂ω ∂y = 0, y = 0, y = H y , 0 ≤ x ≤ H x , 0 ≤ z ≤ H z , где H x = L/ d -длина, H z -высота, H y -ширина
безразмерной расчетной области, R = 0.5 -радиус круглого отверстия; на выходной границе задается условие неотражения [16].

Метод решения

В работе [START_REF] Бекетаева | Ж. // Вычислительные технологии[END_REF][START_REF] Bruel | [END_REF] изложена ENO-схема на основе идеи метода Годунова и продемонстрирована применимость схемы к решению задачи сверхзвукового течения многокомпонентного газа в плоском канале с вдувом перпендикулярных струй. В работе [START_REF] Бекетаева | // Прикладная механика и техническая физика[END_REF] разработанный алгоритм обобщается на трехмерный случай, и показана возможность проведения достоверных численных расчетов сложных сверхзвуковых пространственных течений. В соответствии с этим предварительно в пограничном слое вблизи стенки и на уровне струи для более точного учета течения, вводится сгущение сетки с помощью преобразований [START_REF] Bruel | [END_REF]:

ξ = ξ(x), η = η(z ), ζ = ζ (y). ( 2
)
При этом система уравнений (1) в обобщенных координатах запишется в виде

∂ Ũ ∂t + ∂ Ẽ ∂ξ + ∂ F ∂η + ∂ G ∂ζ = ∂ Ẽv2 ∂ξ + ∂ Ẽvm ∂ξ + ∂ Fv2 ∂η + ∂ Fvm ∂η + ∂ Gv2 ∂ζ + ∂ Gvm ∂ζ . ( 3
) где Ũ = 1 J U, Ẽ = ξ x J E, F = η z J F, Ẽv2 = ξ x J E v2 , Ẽvm = ξ x J E vm , Fv2 = η z J F v2 , Fvm = η z J F vm , Gv2 = ζ y J G v2 , Gvm = ζ y J G vm , J = ∂(ξ, η, ζ ) ∂(x, z , y)
-якобиан преобразования, Ẽvm , Ẽv2 -диффузионные члены, содержащие смешанные и вторые производные. В соответствии с принципом построения ENO-схемы, исходная система уравнений формально представляется следующим образом: [START_REF] Бекетаева | Ж. // Вычислительные технологии[END_REF][START_REF] Bruel | [END_REF]. После факторизации одношаговой конечноразностной схемы для интегрирования по времени уравнения (4) получается следующее равенство:

∂ Ũ ∂t + ( Â+ + Â-) ∂E m ∂ξ + ( B+ + B-) ∂F m ∂η +( Q+ + Q-) ∂G m ∂ζ - ∂( Ẽv2 + Ẽvm ) ∂ξ + ∂( Fv2 + Fvm ) ∂η - ∂( Gv2 + Gvm ) ∂ζ = 0. (4) Здесь A = ∂E ∂U , B = ∂F ∂U , Q = ∂G ∂U -матрицы Якоби, ± = R ˆ ξ R -1 = R 1 ± sign( ξ ) 2 R -1 . B± = T ˆ η T -1 = T 1 ± sign( η ) 2 T -1 , Q± = S ˆ ζ S -1 = S 1 ± sign( ζ ) 2 S -1 , E m = Ẽ + E ξ + D ξ , F m = F + E η + D η , G m = G+E ζ +D ζ - модифицированные потоки в узловых точках (i, j, k), со- стоящие из исходных конвективных векторов ( Ẽ, F, G) и добавочных членов высокого порядка точности (E ξ , D ξ , E η , D η , E ζ , D ζ ), подробно описанных в
I + t ( Â+ + Â-) n ∂ ∂ξ A n ξ • - ∂ ∂ξ μξ ∂ ∂ξ 1 ŨI • × I + t ( B+ + B-) n ∂ ∂η B n η • - ∂ ∂η μη ∂ ∂η 1 ŨI • × I + t ( Q+ + Q-) n ∂ ∂ζ Q n ζ •- ∂ ∂ζ μζ ∂ ∂ζ 1 ŨI • Ũn+1 = Ũn + t ∂ Ẽn v22 ∂ξ + ∂ Fn v22 ∂η + ∂ Gn v22 ∂ζ + ∂ ∂ξ (2 Ẽn vm -Ẽn-1 vm ) + ∂ ∂η (2 Fn vm -Fn-1 vm ) + ∂ ∂ζ (2 Gn vm -Gn-1 vm ) -t ( Â+ + Â-) ∂ ∂ξ (E ξ + D ξ ) + ( B+ + B-) ∂ ∂η (E η + D η ) + ( Q+ + Q-) ∂ ∂ζ (E ζ + D ζ ) n . ( 5
) где A ξ = ξ x A, B η = η z B, Q ζ = ζ y Q причем Â+ + Â-= I, I -единичная матрица, μξ = µξ 2 x Re J , μη = µη 2 z Re J , μζ = µζ 2 y Re J .
Журнал технической физики, 2019, том 89, вып. 10 Аппроксимация производных в конвективных членах, членах содержащих добавочные вектора высокого поряда, а также в членах, содержащих вторые и смешанные производные, производится согласно [START_REF] Бекетаева | Ж. // Вычислительные технологии[END_REF][START_REF] Bruel | [END_REF]. Решение системы (5) осуществляется методом расщепления относительно вектора Ũ матричной прогонкой. На рис. 4, a в плоскости yz (сечение x = 17) показана известная из расчетов при умеренных нерасчетностях [START_REF] Бекетаева | // Прикладная механика и техническая физика[END_REF] пара вихрей [START_REF] Khali | Mixing flow characteristics for a transverse sonic jet injecting into a supersonic crossflow[END_REF], которая генерируется вихрем 5 в результате его бокового перетекания. На рисунке слева видно, что величина завихренности ω = | rot(u)| в области присутствия вихря 8 имеет большие значения. Как следует из вычислительных экспериментов, для рассматриваемого случая эта пара вихрей обладает, так же как и для умеренных n = 50 [START_REF] Бекетаева | // Прикладная механика и техническая физика[END_REF], достаточно большой энергией, что позволяет им не только перемещаться далеко вниз по потоку, но и увеличиваться в размерах. Из численных результатов следует, что максимальные размеры бокового вихря здесь наблюдаются в сечении x = 22.85. Аналогично работе [START_REF] Бекетаева | // Прикладная механика и техническая физика[END_REF] увеличение размеров этих вихрей, по-видимому, обусловлено их локализацией в пристенной области, а рост пограничного слоя проводит к тому, что эти вихри также растут. Из численных расчетов следует, что, начиная с этого сечения, наблюдается уменьшение интенсивности боковых вихрей 8, а также конец бочкообразной структуры в струе. По мере продвижения вихря 8 вниз по потоку заметна потеря его интенсивности в поперечном направлении, поэтому он не наблюдается в сечениях от x = 23.24 до x = 26.90.

Анализ результатов

Ниже на рис. 4, b в плоскости yz (x = 26.30) приведена вихревая пара 9, которая формируется у основания диска Маха в результате взаимодействия струи и восходящего потока под струей. Здесь наблюдается вихревой след 6 вблизи стенки в области низкого давления. Величина завихренности на графике в области вихрей 9 и 6 также имеет большие значения. Наличие этих вихревых пар и механизмов их формирования изучено в работе [START_REF] Бекетаева | // Прикладная механика и техническая физика[END_REF] для умеренных параметров нерасчетности в диапазоне 10 ≤ n ≤ 50, а также в работе [START_REF] Viti | [END_REF] для n = 282. Таким образом, вышеописанная вихревая система является достаточно устойчивой. Распределение величины завихренности показывает присутствие подковообразного вихря 7 в зоне умеренно больших значений завихренности, локализованных в области, примыкающей к стенке.

Картина, где вновь интенсифицированный вихрь 8 вращается в паре с вихрем 9, представлена на рис.

Заключение

Проведен детально-сравнительный анализ механизмов образования вихревых структур в результате взаимодействия набегающего потока с вдуваемой струей при умеренных и больших параметрах нерасчетности 50 ≤ n ≤ 300. В результате выявлено, что, начиная с n = 50, происходит перестройка ориентации центров вращения вихрей 4 и 5 в направлении поперечного потока, следствием которой является затруднительная их визуализация в сечениях, поперечных натекающему Журнал технической физики, 2019, том 89, вып. 10 потоку, а также разделение подковообразного вихря 4 на две независимые вихревые структуры по мере удаления его от линии симметрии. Подтверждено наличие роговых вихрей вблизи стенки в области за струей в диапазоне нерасчетностей 10 ≤ n ≤ 300, возникающих в результате взаимодействия истекающей струи с внешним потоком, приводящего к закрутке встречных потоков за струей в форме двух рогообразных структур. В ходе построения обобщающей зависимости угла наклона головного скачка уплотнения от параметра нерасчетности показан умеренно линейный рост угла при больших значениях нерасчетности.

Финансирование работы

Работа выполнена при поддержке Министерства образования и науки Республики Казахстан в рамках грантового финансирования научных исследований по теме: " Численное моделирование пространственных турбулентных сжамиемых течений с вдувом струей и твердых частиц", 2018-2020 гг., ИРН проекта AP05131555.

Конфликт интересов

Авторы заявляют, что у них нет конфликта интересов.

Список литературы

Рис. 1 .

 1 Схема течения. Пояснения в тексте. Константы принимают следующие значения: σ k = 0.5, σ ω = 0.5, β * = 0.09, β = 0.075, γ * = 5/9.

Исходная система ( 1 )

 1 записана в безразмерной форме. В качестве определяющих параметров приняты размерные параметры на входе û∞ , ρ∞ , T∞ , давление и полная энергия отнесены к значению ρ∞ û2 ∞ , характерным размером длины является диаметр круглого отверстия струи d. Pr -число Прандля, Re = û∞ ρ∞ L/µчисло Рейнольдса, где L -длина канала, являющаяся определяющим размером для Re.

Рис. 2 .Рис. 3 .

 23 a поле вектора скорости и коэффициент давления на стенке CP ; b -проекции линий тока на плоскость, нормальную к оси x; c -линии изобар; d -величина завихренности, при Re = 1.87 • 10 6 , Pr = 0.72, n = 300, M 0 = 1, M∞ = 4, в плоскости симметрии сечение y = 16. Журнал технической физики, 2019, том 89, вып. 10 А.О. Бекетаева, P. Bruel, А.Ж. Найманова X Поле вектора скорости и коэффициент давления на стенке CP (слева), проекции линий тока на плоскость, нормальную к оси y (справа) при Re = 1.87 • 10 6 , Pr = 0.72, n = 300, M 0 = 1, M∞ = 4 в плоскости z x: a -сечение y = 18.41, by = 20.97, cy = 22.98.

4 .

 4 Исследуется обтекание звуковой струи совершенного газа (M 0 = 1) с диаметром отверстия d = 0.476 cm и T0 = 261 K с нерасчетностью в диапазоне 50 ≤ n ≤ 300, сверхзвуковым потоком с параметрами: Pr = 0.72, Re = 1.87 • 10 6 , M ∞ = 4, T∞ = 70.3 K в канале длиной 27.69 cm, высотой 11.43 cm и шириной 15.24 cm соответственно. Безразмерные параметры вычислительной области следующие: H x = 58, H y = 32, H z = 24, где расстояние от входной границы до центра вдуваемой струи составляет x 0 = 16, y 0 = 16 калибров. При этом высота пограничного слоя на входе δ 1 = 1.65, где для этой высоты используется 22-26 узлов расчетной сетки. Диапазон рассматриваемого параметра нерасчетности варьировал в пределах 50 ≤ n ≤ 300. В работе [13] на примере решения пространственной задачи поперечного вдува звуковой струи с умеренно большими параметрами нерасчетности приведен подробный анализ сходимости численного решения по предлагаемой методике. В соответствии с этим в настоящей работе расчетная сетка была выбрана размером 251 × 121 × 151. Ниже на рис. 2-6, продемонстрированы результаты численного моделирования с параметром нерасчетности n = 300. Представленная на рис. 2, a картина поля вектора скорости показывает, что здесь аналогично умеренно большим нерасчетностям формируются два противоположно вращающихся вихря 4 и 5 перед струей, возникающих в результате отрыва набегающего потока, вызванного λ-образной системой скачков уплотнения. Картина проекции линий тока на плоскость, нормальную к оси y в плоскости симметрии (y = 16), хорошо демонстрирует эти вихри (рис. 2, b). На рис. 2, c представлены изобары Журнал технической физики, 2019, том 89, вып. 10 Детально-сравнительный анализ взаимодействия сверхзвукового потока с поперечной газовой струей... Поле вектора скорости и проекции линий тока на плоскость, нормальную к оси x (слева), величина завихренности (справа) при Re = 1.87 • 10 6 , Pr = 0.72, n = 300, M 0 = 1, M∞ = 4 в плоскости z y: a -сечение x = 17, bx = 26.30, cx = 28.18, dx = 42.32. λ-образной системы ударных волн (головной, косой и замыкающий скачки уплотнения). На рис. 2, a также приведен график коэффициента давления на стенке C P (сплошная линия), из которого следует, что по мере приближения к области вдува давление увеличивается, затем его рост замедляется таким образом, что при этом образуется небольшое плато (область I). Далее видно, что давление уменьшается (область II), при этом из поля вектора скорости следует, что область давления II соответствует ядру

6 .

 6 4, c в сечении x = 28.18. Следует отметить, что для параметра нерасчетности n = 282 авторы [9] получили эту структуру с двумя вихрями в сечении x = 31, тогда как численные результаты для этой же нерасчетности выявили это в сечении x = 26.90. Наблюдается качественное совпадение полученных данных, при этом количественное расхождение можно объяснить различной постановкой входных граничных условий для x и ω в Журнал технической физики, 2019, том 89, вып. 10 Проекции линий тока на плоскость, нормальную к оси z (слева), величина завихренности (справа) при Re = 1.87 • 10 6 , Pr = 0.72, n = 300, M 0 = 1, M∞ = 4 в плоскости xy: a -сечение z = 0.1235, bz = 0.2648, cz = 0.5958. модели турбулентности. Здесь же на графике заметно увеличение размеров вихревой пары 6. Далее начиная с сечения x = 30.05, вихрь 9 увлекается вихрем 8. Об этом свидетельствует тот факт, что направление вращения сформированного вихря совпадает с вращением 8 и сохраняется вплоть до конца расчетной области. На рис. 4, d в плоскости z y (x = 42.32) подтверждено присутствие вихревой пары 10, впервые численно выяв-ленной для n = 282 в [9], а также в [13] для умеренных значений нерасчетности (n = 50). Как известно, этот вихрь формируется в результате взаимодействия струи и высокоскоростного натекающего потока, проходящего над вдуваемой струей, и возникает при нерасчетностях выше 20. Этот же рисунок демонстрирует картину возникновения вихря 11, образованного в результате прилипания струи к плоскости пластины и подробно описанного

  τ x y , τ y y , τ y z , uτ y y + vτ y y + wτ y zq y ,

	1 Re	(µ l + σ k µ t )		∂k ∂y	,	1 Re	(µ l + σ ω µ)	∂ω ∂y	T	.
	Тензоры напряжения и потоки тепла выражаются
	в виде
	τ x x =	2 3	µ Re	2	∂u ∂x		-	∂w ∂z	-	∂v ∂y	,
	τ z z =	2 3	µ Re	2	∂w ∂z	-	∂u ∂x	-	∂v ∂y	,
	τ y y =	2 3	µ Re	2	∂v ∂y		-	∂u ∂x	-	∂w ∂z	,
		τ x z = τ z x =	µ Re		∂u ∂z	+	∂w ∂x	,
		τ x y = τ y x =		µ Re		∂u ∂y	+	∂v ∂x	,
		τ y z = τ z y =	µ Re		∂w ∂x	+	∂v ∂z	,

  По мере удаления от линии симметрии вихри 4 и 5 существенно перестраиваются. Начиная с сечения y = 18.74, просматривается только один вихрь 4. Затем его ядро разделяется на два (рис. 3, b, в плоскости xz , y = 20.97). Здесь, как следует из рисунка, максимальное значение C P находится не в зоне разделения двух вихрей, что позволяет предполагать о раздвоении вихря 4. Разделенный вихрь быстро распадается, о чем свидетельствует рис. 3, c, в плоскости xz (y = 22.98), где присутствует только один вихрь. Далее начиная с сечения y = 24.25, вихревая структура 4 уже не просматривается, по-видимому, здесь также происходит перестройка ориентации центра вращения и его вытягивание в направлении поперечного потока.

													А.О. Бекетаева, P. Bruel, А.Ж. Найманова
		80
		60
	deg	40
		20
		0	0	30	60	90	120	150	180	210	240	270	300
								n
	Рис. 5. Влияние параметра нерасчетности на угол наклона
	головного скачка уплотнения.
	подковообразного вихря 4, а его минимальное значение
	приходится на зону центра вращения. За областью
	понижения следует резкое возрастание коэффициента
	давления (область III), где он достигает своего макси-
	мального значения (область IV). Вычислительные расче-
	ты показывают, что максимум давления является зоной
	разделения вращающихся вихрей 4 и 5 (рис. 2, a). Со-
	ответствие максимума давления зоне разделения вихрей
	обусловлено тем, что поток воздуха, проходя участки
	λ-образной системы скачков уплотнения, в частности
	замыкающего скачка, за которым давление достигает
	наибольшего значения, поворачивает вниз к стенке и
	проникает в зону отрывного течения, затем растекается
	в противоположных направлениях с образованием двух
	вихрей, причем вихрь 4 движется по часовой стрелке,
	а 5 -против. На графиках распределения изобар
	(рис. 2, b) видно, что максимум давления находится за
	нижней частью замыкающего скачка уплотнения (значе-
	ние давления 1.5). График 2, d демонстрирует распреде-
	ления величины завихренности ω = | rot(u)|, где отчет-
	ливо видна бочкообразная структура, возникающая при
	вдуве недорасширенной струи в сверхзвуковой поток,
	так же на рисунке можно наблюдать высокие значения
	завихренности в области присутствия вышеописанных
	вихрей.
	Из графика рис. 3, a слева (поле вектора скорости и
	коэффициент давления на стенке C P в плоскости xz ,
	y = 18.41) можно наблюдать боковое растекание вих-
	рей 4 и 5. Из рисунка видно, что оба вихря разделе-
	ны максимальным значением давления на графике C P
	(область IV). Из рисунка следует нечеткая вихревая
	структура 5 (рис. 3, a справа -проекции линий тока на
	плоскость, нормальную к оси X), заметно лишь ее уве-
	личение восходящим потоком. По всей видимости здесь
	происходит перестройка ориентации центра вращения
	вихря 5.

  Распределение коэффициента давления CP на стенке на оси симметрии при Re = 1.87 • 10 6 , Pr = 0.72, n = 300, M 0 = 1, M∞ = 4: сплошная линия -численные результаты, -эксперимент[START_REF] Viti | [END_REF], ----результаты работы[START_REF] Viti | [END_REF].расхождение расчетных и экспериментальных данных в диапазоне -4 < x < -1, по всей видимости, обусловлено использованием k-ω модели турбулентности, которое, возможно, повлияло на область стыковки вдуваемой струи и основного потока вблизи стенки, что свидетельствует о том, что необходима дальнейшая модификация данной модели. Этот факт также подтверждает работа[START_REF] Волков | // Прикладная механика и техническая физика[END_REF], которая посвящена струйному взаимодействию. В настоящей работе изучена применимость различных моделей турбулентности, а именно модели Спаларта-Альмараса, k-ε, k-ω и SST-моделей при решении подобных задач на основе сравнения вычислительных и экспериментальных данных. Здесь показано, что использование различных моделей турбулентности приводит к различным численным результатам по распределению давления на стенке перед вдуваемой струей. В частности, расчет с k-ω моделью турбулентности приводит к занижению результатов численного моделирования по сравнению с экспериментом также в области, близкой к вдуваемой струе.

						А.О. Бекетаева, P. Bruel, А.Ж. Найманова
	в [13] и, как следует из рисунка, направление его вращения совпадает с вращением следа 6. Численные		0.5
	эксперименты показывают, что в дальнейшем вихри 11 и 6 сливаются, и вниз по потоку наблюдается только		0.4
	один вихревой след. Из картины распределения завих-
	ренности видно, что подковообразный вихрь 7 в области		0.3
	умеренных значений отходит от линии симметрии на
	большое расстояние. Следует отметить, что локализация вихря 10 по	C p	0.2
	высоте напрямую зависит от интенсивности вдува струи
	и ударно-волновой структуры. Так, при умеренных зна-		0.1
	чениях нерасчетности (см. работу [13] при n = 50),
	данный вихрь формируется на высоте от 4 до 5 калиб-ров. В рассматриваемом случае локализация его центра		0
	вращения оказалась на уровне z = 8.7, а при нерасчет-
	ности n = 282 -на высоте z = 8, 5. Так как вихрь (10) ком уплотнения, можно сделать вывод, что его высота расположен между вдуваемой струей и головным скач-		-12	-8	-4	X 0	4	8	12
	напрямую зависит от угла наклона головной ударной	Рис. 7.
	волны. На рис. 5 показан график влияния параметра
	нерасчетности на головной скачок уплотнения, где ось
	абсцисс соответствует нерасчетности, ось ординат -уг-
	лу наклона головного скачка, вычисленного по формуле
	ang = arctg max	w u	.
	Здесь видно, что угол наклона головного скачка уплот-
	нения возрастает с увеличением отношения давления
	в струе и потоке. Как видно из графика, угол наклона
	при малых и умеренных значениях нерасчетности резко
	возрастает, однако, начиная с n = 150, значения стано-
	вятся практически постоянными, и высота локализации
	вихрей 10 будет меняться незначительно, т. е. начиная с
	n = 150, высота этих вихрей практически постоянна.
	На рис. 6, a в плоскости xy (z = 0.1235) можно отчет-
	ливо увидеть вихри 12, которые сформированы вблизи
	стенки в области за струей на небольшом расстоянии
	от отверстия (в сечении x = 19.5). Данные вихри пред-
	ставлены в работе [10] и носят название horn (роговые
	вихри). Физический механизм их возникновения состоит
	в том, что струя, истекая непосредственно из отверстия,
	взаимодействует с внешним потоком, который, обтекая
	струю, стремится замкнуться на линии симметрии, что
	приводит к закрутке встречных потоков в форме двух
	рогообразных структур. Численные результаты показы-
	вают, что этот вихрь формируется и при параметре
	нерасчетности n = 5, что свидетельствует об устойчи-
	вости данной вихревой системы. Численные результаты
	показывают, что в сечении z = 0.2648 (рис. 6, b) вихри
	замыкаются на линии симметрии, а уже при z = 0.5958
	(рис. 6, c) роговые вихри не просматриваются.
	На рис. 7 представлено сравнение численного резуль-
	тата коэффициента давления C P на стенке в плоскости
	симметрии с данными работы [9]. Как видно из рисунка,
	получено удовлетворительное согласование численных
	результатов и экспериментальных данных. Небольшое

Журнал технической физики, 2019, том 89, вып. 10

Журнал технической физики, 2019, том 89, вып. 10

* Журнал технической физики, 2019, том 89, вып. 10


Mohamed Elshamy

Alhassan Sayed

Marie-Minerve Louërat

Amine Rhouni

Hassan Aboushady

Haralampos-G Stratigopoulos

Securing Programmable Analog ICs Against Piracy

In this paper, we demonstrate a security approach for the class of highly-programmable analog Integrated Circuits (ICs) that can be used as a countermeasure for unauthorized chip use and piracy. The approach relies on functionality locking, i.e. a lock mechanism is introduced into the design such that unless the correct key is provided the functionality breaks. We show that for highly-programmable analog ICs the programmable fabric can naturally be used as the lock mechanism. We demonstrate the approach on a multi-standard RF receiver with configuration settings of 64-bit words.

I. INTRODUCTION

Hardware security and trust is a topic that has attracted a lot of interest in recent years. There are various threats, including IC/IP piracy, hardware Trojans, side-channel attacks, and fault injection attacks, which can have serious implications on the technological value chain (e.g. CAD tool providers, IC/IP providers, original equipment manufacturers, and users), on governments, and on the society as a whole [START_REF] Rostami | A primer on hardware security: Models, methods, and metrics[END_REF].

While hardware security and trust aspects have been extensively studied for digital circuits, the space of vulnerabilities and solutions for analog circuits is largely unexplored and little understood as of today [START_REF]Trusted analog/mixed-signal/RF ICs: A survey and a perspective[END_REF], [START_REF] Alam | Challenges and Opportunities in Analog and Mixed Signal (AMS) Integrated Circuit (IC) Security[END_REF].

In this paper, we address the problem of analog IC/IP piracy, which includes reverse engineering and counterfeiting. Reverse engineering refers to the derivation of IC/IP proprietary information, i.e. architecture, netlist, layout, etc. It aims at reducing the attacker's technological disadvantage against the "author" of the IC/IP, gathering necessary information for producing a similar or identical IC/IP, e.g. a counterfeit, or locating the root-of-trust part of the IC/IP to steal secret information, such as cipher keys. Nowadays, there exist equipment and software tools to successfully reverse-engineer any unprotected IC/IP [START_REF] Torrance | The state-of-the-art in semiconductor reverse engineering[END_REF]. Counterfeiting includes cloning, recycling, overproducing, and remarking [START_REF] Guin | Counterfeit integrated circuits: A rising threat in the global semiconductor supply chain[END_REF]. A cloned counterfeit is an IC/IP that is illegally cloned and sold as original. Cloning can be performed by an untrusted foundry or an adversary via reverse-engineering. A recycled counterfeit is a used and possibly aged IC that is illegally resold as new. Overproduced ICs are ICs that are produced by an untrusted foundry beyond the number agreed in the contract with the IC design house and are illegitimately sold in after market. Remarked ICs are failing ICs that are remarked by an untrusted test facility as passing ICs and are sold with false and forged documentation.

In particular, in this paper, we propose a security approach based on locking for the class of highly-programmable analog ICs. The approach can be used not only for prevention of unauthorized chip use, but, in addition, it can offer a strong countermeasure against reverse engineering and counterfeiting. [START_REF] Hoe | Towards secure analog designs: A secure sense amplifier using memristors[END_REF]; (b) obfuscating biasing transistors [START_REF] Rao | Protecting analog circuits with parameter biasing obfuscation[END_REF]; (c) locking of current mirrors [START_REF] Wang | Thwarting analog IC piracy via combinational locking[END_REF]; (d) locking mixed-signal circuits via logic locking of their digital section [START_REF] Leonhard | Mixlock: Securing mixed-signal circuits via logic locking[END_REF]; (e) logic locking of the digital optimizer in the calibration feedback loop [START_REF] Jayasankaran | Towards provably-secure analog and mixedsignal locking against overproduction[END_REF]; (f) locking through neural network-based biasing [START_REF] Volanis | Analog performance locking through neural networkbased biasing[END_REF].

We argue that the tuning knobs within the programmable analog IC can naturally serve as a locking mechanism. Specifically, the programming bits controlling the tuning knobs serve as key-bits and each configuration setting, e.g. programming bits that configure the IC in a specific operation mode demanded by the application, is treated as a secret key. Naturally, when invalid programming bits are provided the functionality of the circuit breaks, that is, its functionality is "locked". We discuss the practical implementation of this security approach, its benefits compared to existing locking techniques, and its resilience against foreseen attacks. We demonstrate it on a highly-programmable, multi-standard RF receiver with configuration settings of 64-bit words.

The rest of the paper is structured as follows. In Section II, we discuss previous work on analog circuit locking. In Section III, we provide a general overview of programmability embedded into analog ICs. In Section IV, we present the proposed technique for securing programmable analog ICs. In Section V, we present our case study. In Section VI, we demonstrate the efficiency of locking and the achieved security level. Section VII concludes the paper.

II. PREVIOUS WORK ON ANALOG CIRCUIT LOCKING

Techniques for locking of analog circuits are proposed in [START_REF] Hoe | Towards secure analog designs: A secure sense amplifier using memristors[END_REF]- [START_REF] Volanis | Analog performance locking through neural networkbased biasing[END_REF] and are illustrated in Fig. 1. In [START_REF] Hoe | Towards secure analog designs: A secure sense amplifier using memristors[END_REF], a locking mechanism based on an architecture that comprises memristor crossbars is used to lock the body biasing of the transistor input pair in a sense amplifier. In [START_REF] Rao | Protecting analog circuits with parameter biasing obfuscation[END_REF], it is proposed to replace transistors within the biasing circuit with parallel-connected transistors whose gates are controlled by key-bits. The keybits set an aggregate width equal to the width of the original transistor. In [START_REF] Wang | Thwarting analog IC piracy via combinational locking[END_REF], it is shown how to redesign the current mirrors providing the biasing so as to insert key-bits. In [START_REF] Leonhard | Mixlock: Securing mixed-signal circuits via logic locking[END_REF], it is proposed to lock a mixed-signal circuit via logic locking of its digital section. This technique was demonstrated in an audio application in [START_REF] Leonhard | Mixed-signal hardwre security using mixlock: Demonstration in an audio application[END_REF] allowing to listen to the effect of locking on the audio quality. In [START_REF] Jayasankaran | Towards provably-secure analog and mixedsignal locking against overproduction[END_REF], it is proposed to lock via logic locking the digital optimizer into the calibration feedback loop, such that the wrong tuning settings are generated unless the valid key is applied. In [START_REF] Volanis | Analog performance locking through neural networkbased biasing[END_REF], it is proposed to add on-chip a neural network that is trained to map the secret analog key, which is in the form of analog DC voltages presented as inputs to the neural network, to the correct biases.

The locking approaches in [START_REF] Hoe | Towards secure analog designs: A secure sense amplifier using memristors[END_REF]- [START_REF] Wang | Thwarting analog IC piracy via combinational locking[END_REF], [START_REF] Volanis | Analog performance locking through neural networkbased biasing[END_REF] act on the biasing of circuit. This makes them vulnerable to removal attacks since the attacker does not have to recover the key; it suffices to recover the biases, which are typically small in number, and thereafter replace the locked blocks with "fresh ones" that provide the correct biases.

In contrast to biases that are fixed for every fabricated chip, tuning knobs are set per fabricated chip to compensate for process variations, thus to break the technique in [START_REF] Jayasankaran | Towards provably-secure analog and mixedsignal locking against overproduction[END_REF] the attacker will need to actually search for the secret key. The same holds for the technique in [START_REF] Leonhard | Mixlock: Securing mixed-signal circuits via logic locking[END_REF] which does not lock the biasing or tuning knobs, but locks directly the functionality. Still, [START_REF] Leonhard | Mixlock: Securing mixed-signal circuits via logic locking[END_REF], [START_REF] Jayasankaran | Towards provably-secure analog and mixedsignal locking against overproduction[END_REF] are vulnerable to removal attacks as the attacker may replace the locked digital optimizer in [START_REF] Jayasankaran | Towards provably-secure analog and mixedsignal locking against overproduction[END_REF] or the locked digital section in [START_REF] Leonhard | Mixlock: Securing mixed-signal circuits via logic locking[END_REF] with "fresh" unlocked designs, although arguably this is a much more difficult task than redesigning the biasing circuitry.

III. PROGRAMMABLE ANALOG ICS

Analog circuits are often made programmable (or configurable) with the aim to: (a) Compensate for process variations so as to increase yield; (b) Compensate for inherent nonidealities so as to achieve the desired performance tradeoff; (c) Configure the circuit into different operation modes demanded by the application; (d) Adapt the performances to changes in the environment.

The calibration mechanism that enables programmability (or configuration) consists of tuning knobs judiciously inserted into the design and a calibration algorithm that is driven by performance indicators and returns the programming bits (or configuration setting). Often the same calibration mechanism is used to achieve simultaneously multiple of the above objectives. Typically process variations are taken into consideration during calibration, thus the configuration settings end up being unique for each chip.

The calibration algorithm can either run on-chip in hardware pointing to autonomous self-calibration or can run off-chip in software during the post-manufacturing testing phase, relying on external Automated Test Equipment (ATE). For objectives (c)-(d), operation modes and adaptation levels are pre-specified based on the anticipated range of applications and environmental conditions, resulting in multiple pre-specified configuration settings that are pre-loaded into an on-chip look-up table (LUT). Instead, for objectives (a)-(b), the calibration algorithm returns a single best possible configuration setting and the LUT has essentially one line.

The above discussion remains quite general; in fact, the calibration mechanism varies from one circuit class to another, and programmability may vary from a few bits for calibrating single blocks [START_REF] Andraud | One-shot non-intrusive calibration against process variations for analog/RF circuits[END_REF] to tens of bits for calibrating complete systems [START_REF] Li | Reconfigurable All-Band RF CMOS Transceiver for GPS/GLONASS/Galileo/Beidou With Digitally Assisted Calibration[END_REF].

IV. SECURING HIGHLY-PROGRAMMABLE ANALOG ICS

A. Locking via the programmability fabric

In this work, we argue that for highly-programmable circuits it is not required to insert additional circuitry on-chip in order to introduce key-bits. Instead, we can take advantage of the embedded programmable fabric so as to naturally perform the locking operation, as shown in Fig. 2. Specifically, the configuration settings can be treated as secret keys or, equivalently, the programming bits can be treated as secret key-bits. Similarly, the calibration algorithm that produces the configurations settings is kept secret and is not shared with untrusted parties. Using invalid programming bits will result in complete loss of functionality or in significant performance degradation, that is, one or more performances will lie far outside their allowable specification range.

Exploiting the embedded programmability for performing the locking operation presents significant advantages. Unlike known approaches for locking analog ICs which all modify to some degree the design [START_REF] Hoe | Towards secure analog designs: A secure sense amplifier using memristors[END_REF]- [START_REF] Volanis | Analog performance locking through neural networkbased biasing[END_REF], with the proposed approach the design is left completely intact. Therefore, there is no need for redesign, no extra design iterations, and, most importantly, no performance degradation. In addition, the proposed approach does not increase the power or area of the analog IC itself. The power and area overheads are only due to the key management scheme, which can be shared for enabling security for all other blocks on the same die, i.e. in the context of a System-on-Chip (SoC).

The possible secret key management schemes are the same ones used by logic locking techniques for digital ICs [START_REF] Roy | Ending piracy of integrated circuits[END_REF]. One option is to store the LUT with the configuration settings into a tamper-proof memory, as shown in Fig. 3(a). A second option, illustrated in Fig. 3(b), makes use of a Physical Unclonable Function (PUF) [START_REF] Herder | Physical unclonable functions and applications: A tutorial[END_REF]. The PUF needs to take at least as many challenges as the total number of configuration settings. The PUF generates a number of secret identification keys that equals the number of configuration settings. The user is given a number of keys that also equals the number of configuration settings, such that when these keys are XORed with the identification keys the correct configuration settings are produced. In both schemes, in normal operation mode the circuit commands dynamically the memories to load the corresponding programming bits.

B. Attack scenarios and resilience analysis

We assume that an attacker has full capabilities, i.e., has the netlist and access to working oracle chips.

1) Attacks in digital domain: Known attacks in digital domain, such as the lethal SAT attack [START_REF] Subramanyan | Evaluating the security of logic encryption algorithms[END_REF], are not applicable.

2) Removal attacks: Removal attacks, which is the main limitation in [START_REF] Hoe | Towards secure analog designs: A secure sense amplifier using memristors[END_REF]- [START_REF] Volanis | Analog performance locking through neural networkbased biasing[END_REF] as described in Section II, are not applicable as there is no added circuitry on-chip to facilitate the key insertion; the key directly applies to existing tuning knobs into the design.

3) Brute-force and multi-objective optimization attacks: The most trivial attack is the brute-force attack which consists in applying random combinations of programming bits until the one that unlocks the circuit is found.

The multi-objective optimization attack consists in applying an iterative algorithm that searches for a configuration setting that simultaneously optimizes the performances such that they all satisfy their specifications. This attack is difficult to put in place since typically only a small subset of programming bits shows a smooth monotonic relationship with a given performance, and this requires that the rest of the programming bits are already correctly set.

Performing the brute-force and multi-objective optimization attacks by simulation is impractical due to very long analog simulation times. To perform these attacks much faster in hardware, the attacker needs first to re-fab the circuit so as to gain direct access to programming bits. In this case, if the programming bits are unique for each chip, then these attacks become meaningful only if the resultant key-bit combination can be used to set a good starting point for launching a gradient search for quickly calibrating any chip.

A question rises whether the design can be divided in subblocks, tracing key bits to sub-blocks, and enabling smaller brute-force and multi-objective optimization attacks at subblock level. This is typically not possible due to the internal feedback loops that involve multiple sub-blocks each. Moreover, this would require re-fabbing so as to provide intermediate taps to access directly sub-blocks for measurement, which is only possible for low-frequency sub-blocks.

It should be noted that unlike logic locking of digital circuits where there is a single valid key, for analog circuits it is likely that a number of key-bit combinations result in a satisfactory performance trade-off, although this number is typically a very small fraction of all key-bit combinations.

Clearly, resilience against these attacks increases with the number of programming bits and with simulation or measurement time per trial.

4) Revealing the calibration algorithm: In the case of an on-chip calibration, an attacker that extracts the netlist will also have at hand the hardware that implements the calibration feedback loop and, thereby, it may be fairly easy to extract the calibration algorithm. In this scenario, we can envision logic locking of the digital section of the calibration feedback loop [START_REF] Leonhard | Mixlock: Securing mixed-signal circuits via logic locking[END_REF], [START_REF] Jayasankaran | Towards provably-secure analog and mixedsignal locking against overproduction[END_REF].

In the case of an off-chip calibration, the attacker may target speculating the calibration algorithm by studying the circuit architecture. However, very often the calibration algorithm is very specific and esoteric to the design, that is, intended for or understood only by the designer. Thus, the attacker must have a very high and specialized expertise and a thorough understanding of the design so as to be able to conceive the underlying calibration algorithm. This is in contrast to the security assumption, in the sense that if the attacker is so knowledgeable, then the attacker may as well design the circuit from scratch instead of pirating an existing one, given also that piracy by an attacker other than the foundry requires a sophisticated and expensive reverse-engineering infrastructure. This is a new type of attack specific to the countermeasure that is proposed. In general, it opens a discussion of securing and obfuscating calibration algorithms when they are considered to be a valuable intellectual property of the design. A metric to quantify the difficulty for reverse-engineering a calibration algorithm will need to devised also.

A problem that needs to be addressed is how to protect the calibration algorithm against an untrusted test facility. For programmable analog circuits, multiple test/calibration iterations are carried out during testing for searching in the space of configuration settings. The search is not random, but is driven by the calibration algorithm. In each step of the algorithm, the next configuration setting is dictated by tests done using the current configuration setting. If the calibration runs entirely on-chip, then it is transparent to the test facility. If the calibration runs off-chip, then the obvious workaround is that the foundry returns the manufactured chips to the design house, where the test/calibration steps are performed in a secured environment. However, this option is practical only in the case of low-volume products. For high-volume products, it is straightforward to adapt the concept of remotely activating the chips using asymmetric cryptography [START_REF] Roy | Ending piracy of integrated circuits[END_REF].

C. Threats being addressed

Assuming resilience against the aforementioned attacks, the proposed approach can be used as countermeasure for reverseengineering, cloning, overproducing, and remarking. It can offer resilience against recycling only if the key management scheme in Fig. 3(b) is used and the unique user keys are loaded every time at power-on. More specifically, a cloned counterfeit that can be produced by reverse-engineering a chip or by the untrusted foundry is good-for-nothing if the adversary does not know how the design can be programmed. Resilience against overproducing can be achieved since the design house can have control over the number of correctly programmed fabricated chips. Remarking can be achieved since the design house can load after unsuccessful calibration wrong configuration settings to render the chip totally malfunctional.

V. CASE STUDY: PROGRAMMABLE MULTI-STANDARD RF RECEIVER

Our case study is a programmable highly-digitized multistandard 65 nm CMOS RF receiver. It is made re-configurable such that it can serve for establishing communication using several standards within the frequency range from 1.5 GHz to 3.0 GHz, including Bluetooth, ZigBee, WiFi 802.11b, etc.

A. Architecture and Programmability

Fig. 4 shows the block-level schematic of the RF receiver. It is composed of a Variable Gain Low Noise Amplifier (VGLNA), a band-pass (BP) RF Σ∆ modulator, a digital down-conversion mixer, and a digital decimation filter. In total, there are 64 programming bits embedded into the analog section and 3 programming bits embedded into the digital section. The calibration is run off-chip and returns the programming bits per standard in the presence of process variations, thus resulting in unique configuration settings per standard and per chip. For the purpose of locking, we consider only the programming bits of the analog section since the calibration of the digital section for a given standard is straightforward.

The block-level schematic of the VGLNA is shown in Fig. 5. It is composed of five gain stages with a resistive feedback. It features a 4-bit configuration word with which the VGLNA can attain 16 different gain levels so as to adapt the sensitivity and dynamic range of the RF receiver to the specifications imposed by the target standard.

The block-level schematic of the BP RF Σ∆ modulator is illustrated in Fig. 6 [START_REF] Ashry | A 4th order 3.6GS/s RF Sigma-Delta ADC with a FoM of 1pJ/bit[END_REF]. It is composed of an input transconductance, Gm in , an LC bandpass loop filter with two capacitor arrays C c and C f for coarse-and fine-tuning, a pre-amplifier, a comparator, a loop delay, a feedback Digital-to-Analog Converter (DAC), and an output buffer. For a target standard, the modulator uses a 60-bit configuration word to tune the center frequency and quality factor of the LC bandpass loop filter in the presence of process variations, as well as to trim the biasing current of the other blocks in order to compensate for process variations and improve the performance trade-off.

B. Calibration Algorithm

The calibration procedure for a target standard is as follows:

1) The comparator is configured as a buffer by deactivating its driving clock.

2) The output buffer is configured to adapt the output of the BP RF Σ∆ modulator to its off-chip load during calibration. This output buffer is removed from the signal path in normal operation mode.

3) The RF input signal is disabled by turning off the input transconductor Gm in . 4) The feedback loop with the DAC and loop delay is turned off. 5) Having deactivated the feedback loop, the LC loop filter is put in oscillation mode by setting its Q-enhancement transconductor, -ve Gm, to its maximum. 6) The capacitor arrays C c and C f of the LC tank are tuned until the output frequency is equal to the desired center frequency. 7) The Q-enhancement transconductor, -ve Gm, is reduced gradually until oscillation vanishes. 8) The feedback loop is restored. 9) The BP RF Σ∆ modulator is put in the operating mode by applying an RF input signal with frequency F 0 . 10) The sampling frequency is set to

F s = 4 • F 0 .
11) The loop delay is set according to F s . 12) The VGLNA is tuned to set the appropriate sensitivity and dynamic range. 13) The input transconductance Gm in , the feedback DAC, the pre-amplifier and the comparator are initialized to their nominal values determined by simulation. 14) An iterative procedure is used to determine the configuration words of these blocks through the improvement of the measured Signal-to-Noise Ratio (SNR) and Spurious Free Dynamic Range (SFDR) of the BP RF Σ∆ modulator.

VI. RESULTS

A. Locking efficiency

For a given center frequency or standard, there is an optimal combination of the 64 programming bits composing a secret key that unlocks functionality. We will consider the maximum center frequency, e.g. 3 GHz, and we will demonstrate the locking efficiency when applying invalid keys. The circuit has several performances, including SNR, dynamic range, SFDR, etc., and locking succeeds when at least one performance violates its specification.

We assume that the attacker has extracted the netlist of the circuit and can simulate it at transistor-level with the ability to monitor internal nodes that shed more light into the operation. We consider first the SNR observed at the output of the BP RF Σ∆ modulator for an input sinusoidal signal with frequency 3 GHz and power -25 dBm. The SNR is computed for an Oversampling Ratio (OSR) of 64 and based on a 8192 point FFT. Fig. 7 shows the SNR across 100 randomly keys and the correct key. As it can be seen, the correct key stands out resulting in an SNR of over 40 dB, while for invalid keys the SNR is less than 30 dB. In fact, for most invalid keys the SNR is below 0 dB, which means that the input signal gets buried under the noise level or there are harmonics within the band-of-interest. However, there are 4 invalid keys that have an SNR higher than 10 dB, and among them one has an SNR of about 30 dB. This key with index 7 in Fig. 7 naturally would attract the attention of the attacker. But, in fact, it turns out to result in a "deceptive" SNR. Looking into the programming bits corresponding to this key, the feedback loop of the BP RF Σ∆ modulator is open and, in addition, the comparator operates as a buffer. In this way, the analog signal passes without being digitized, thus there is no quantization noise being added. Fig. 8 shows the transient output of the BP RF Σ∆ modulator for the correct key and the invalid key with index 7. The correct output is an oversampled bitstream, while the output for the invalid key is an analog waveform showing no analog-to-digital conversion. This analog waveform when it passes directly through the digital section of the RF receiver will show a reduced SNR. This is shown in Fig. 9, where now the SNR measured at the output of the RF receiver is plotted. The SNR for the correct key does not change as expected and for some invalid keys the SNR is further reduced. In short, all invalid keys show an SNR of less than 10 dB, that is, the functionality is significantly corrupted. Fig. 10 shows the Power Spectral Density (PSD) at the output of the BP RF Σ∆ modulator for the correct key and the invalid key with index 7. As it can be seen, for the invalid key there is no noise shaping, which is the main characteristic of the BP RF Σ∆ modulator. Fig. 11 shows the dynamic range of the RF receiver for the correct key and the invalid key with index 7. The input range is divided into three segments, e.g. [-85:-45], [-60:-20], and [-40:0], and for each segment the VGLNA is tuned to have the appropriate gain level and sensitivity. While Figs. 7 and9 show the SNR just for an input power of -25 dBm, Fig. 11 plots the SNR for different input power values with a step of 5 dBm. As it can be seen, the behavior of the locked circuit across the input range is very different as compared to the unlocked circuit.

Fig. 12 shows the SFDR for the correct key and the invalid key with index 7. SFDR is measured by applying a twotone input, where the two tones have the same power and a frequency difference of 10 MHz. SFDR is the difference between the power of the fundamental and the third harmonic. As it can be seen, the locked circuit has a much lower SFDR.

Finally, the same experiment was repeated for other center frequencies and qualitatively the results were identical.

B. Security Analysis

1) Brute-force and multi-objective optimization attacks: Simulation of this device is too time-consuming, thus the attacker will need to re-fab the chip so as to perform the analysis in hardware. For example, for a single key and a 8192 point FFT, it takes about 20 minutes to simulate the SNR at the output of the RF receiver for a given input, 3 hours to simulate the SNR across the input range, and 30 minutes to simulate the SFDR.

Most sub-blocks are included in a feedback loop which makes it impossible to calibrate individual sub-blocks. Also, to calibrate a sub-block, the rest of the sub-blocks need to be conditioned appropriately. Besides, the circuit carries highfrequency signals and re-fabbing a chip with intermediate taps for observing internal signals would result in performance loss.

Resilience against these attacks is also naturally achieved thanks to the large key-width of 64 bits which explodes the search space. It is very unlikely that many key-bit combinations could result in satisfactory performance trade-off. For example, capacitor arrays are binary-weighted, thus for a desired capacitor value there is a unique sub-key.

2) Speculating the calibration algorithm: The calibration procedure is arguably very complex and its steps cannot be easily retraced by conjecture even under the assumption that the attacker has strong analog design expertise. In particular: (a) The circuit needs to be reconfigured appropriately multiple times during calibration; (b) The calibration of many subblocks requires initial programming bits that are dictated by design simulation and are unknown to the attacker; (c) The order with which the different blocks should be calibrated is very specific; (d) The feedback loop prohibits calibration of individual sub-blocks.

VII. CONCLUSION

We propose a locking methodology for highlyprogrammable analog ICs that relies on securing the configuration settings as well as the underlying calibration algorithm. The methodology can be used as a countermeasure for unauthorized chip use and piracy. It is demonstrated on a multi-standard RF receiver with 64 programming bits. The proposed methodology has the advantage that it is strictly non-intrusive to the design and has zero power and area overheads, except for the overheads due to the key management scheme.

Fig. 1 .

 1 Fig. 1. Locking techniques for analog ICs: (a) locking biases based on memristor crossbars [6]; (b) obfuscating biasing transistors [7]; (c) locking of current mirrors[START_REF] Wang | Thwarting analog IC piracy via combinational locking[END_REF]; (d) locking mixed-signal circuits via logic locking of their digital section[START_REF] Leonhard | Mixlock: Securing mixed-signal circuits via logic locking[END_REF]; (e) logic locking of the digital optimizer in the calibration feedback loop[START_REF] Jayasankaran | Towards provably-secure analog and mixedsignal locking against overproduction[END_REF]; (f) locking through neural network-based biasing[START_REF] Volanis | Analog performance locking through neural networkbased biasing[END_REF].

Fig. 2 .

 2 Fig. 2. Locking via the programming fabric.

Fig. 3 .

 3 Fig. 3. Key management schemes.

Fig. 4 .

 4 Fig. 4. Architecture of programmable multi-standard RF receiver.

Fig. 5 .

 5 Fig. 5. Architecture of tunable variable gain LNA.

Fig. 6 .

 6 Fig. 6. Architecture of tunable BP RF Σ∆ modulator.

Fig. 7 .

 7 Fig. 7. SNR for correct key (blue cross) and invalid keys (gray dots and red dot with index 7) computed at the output of the BP RF Σ∆ modulator.

Fig. 8 .

 8 Fig. 8. Transient output of BP RF Σ∆ modulator for the correct key (top) and the invalid key with index 7 in Fig. 7 (bottom).

Fig. 9 .

 9 Fig. 9. SNR for correct key (blue cross) and invalid keys (gray dots and red dot with index 7) computed at the output of the RF receiver.

Fig. 10 .

 10 Fig. 10. PSD at the output of the BP RF Σ∆ modulator for the correct key (blue) and the invalid key with index 7 in Figs. 7 and 9 (red).

Fig. 11 .

 11 Fig. 11. SNR versus input power with different LNA gain settings for the correct key (blue) and the invalid key with index 7 in figures 7 and 9 (red).

Fig. 12 .

 12 Fig. 12. SFDR for the correct key (blue) and the invalid key with index 7 in Figs. 7 and 9 (red).

ACKNOWLEDGMENTS

This work has been carried out in the framework of the ANR STEALTH project with N o ANR-17-CE24-0022-01. It is partially funded by the ANR TOLTECA project with N o ANR-16-CE04-0013-01.


A Inés Fernández

Aran Solé

Jessica Giró-Paloma

Mònica Martínez

Mila Hadjieva
email: hadjieva@phys.bas.bg

Abderrahim Boudenne

Mariaella Constantinescu
email: mariella_const@yahoo.com

Elena Maria Anghel

Marta Malikova
email: marta.malikova@savba.sk

Igor Krupa

A Inés Fernández
email: ana_inesfernandez@ub.edu

Abdel Boudenne

Conchita Peñalosa

Ana Lázaro
email: ana.lazaro@unizar.es@unizar.es

Halime O Paksoy
email: hopaksoy@cu.edu.tr

Kemal Cellat

Jana Vecstaudža

Diana Bajare
email: diana.bajare@rtu.lv

Bostjan Sumiga
email: bostjansumiga@gmail.com

Bojana Boh

Thomas Haussmann
email: thomas.haussmann@ise.fraunhofer.de

Stefan Gschwander

Robert Weber
email: robert.weber@empa.ch

Piotr Furmanski

Maciej Jaworski
email: maciej.jaworski@itc.pw.edu.pl

Luisa F Cabeza
email: lcabeza@diei.udl.cat

Unconventional experimental technologies used for

Keywords: Phase change material (PCM), characterization, morphological and structural characterization, physico-chemical stability, mechanical properties

Introduction

Characterization of PCM and hybrid PCM composites is of great interest for the scientific community, because it is of extreme importance for the deployment of such materials into the market. In a previous paper [START_REF] Cabeza | Unconventional technologies available for phase change materials (PCM) characterization[END_REF] special attention was paid to the characterization of thermophysical properties showing the difficulties of using commercial equipment with composite materials and/or samples size above few milligrams.

Although the thermopysical properties are the main selection criteria for the use of PCM, other important properties such as chemical stability, crystalline structure or mechanical properties are extremely important to predict the performance of the hybrid PCM composites. This paper includes examples for the analysis of morphological and structural characterization of hybrid PCM composites, physico-chemical stability and mechanical properties.

Characterization of hybrid materials has different approaches depending on the scale (nano-and micro-to bulk and macro-scale). The morphological and structural characterization gives information on crystalline characteristics, and how the microstructure may be affected by cycling the hybrid material. PCM crystalline structure is evaluated with cycles and temperature changes. Matrix structure or even the encapsulating coating materials may be also evaluated in terms of volume expansion while melting, and therefore the change in their mechanical properties.

Chemical stability of the hybrid materials is also a key issue. Cycling organic materials may lead to a partial degradation that is reflected then in thermal properties. Related with the chemical stability there is the compatibility of the composite, how it is affecting the PCM the matrix properties i.e. chemical, mechanical or even fire resistance properties.

Morphological and structural characterization of hybrid PCM composites

The morphological and structural characterization of hybrid PCM composites may be carried out using commercial and experimental devices. This paper is focussed in both types, giving examples of the characterization of different PCM and PCM-composites using several experimental procedures.

X-Ray Diffraction

Crystalline structure of organic and inorganic PCM and hybrid PCM composites is studied by X-Ray Diffraction at the Institute of Physical Chemistry Romanian Academy (Romania) [START_REF] Constantinescu | Testing Methods for Characterization of New Materials PCM-Epoxy to be used in Buildings[END_REF].

Changes after cycling may be followed by this technique. Moreover, it is possible to calculate the crystallite size, and with additional data the degree of crystallinity and the time of the fastest crystallization. Figure 1 shows an example where the crystalline structure of a polyethylene glycol-epoxy (PEG-epoxy) composite is evaluated. X-ray diffractograms were collected by means of a Rigaku diffractoctomer type Ultima IV in parallel-beam geometry. The X-ray from a Cu tube operating at 40 kV and 30 mA. Counts were collected from 5 0 to 60 0 with a step size of 0.02 and a speed of 5 0 /min. It was assumed that the main contributor to the peak broadening is the finite size of the PEG chains.

Optical and thermal microscopies

The optical microscopy is used at the Bulgarian Academy of Sciences (Bulgaria) to control crystal formation, size, allocation, and modification before and after thermal cycling in PCM slurries. An example can be seen in Figure 2, where two morphological photos compared are of Clathrate hydrate slurry. "In situ" control during thermal cycling of PCM can be carried out with a thermal microscope (Bulgarian Academy of Sciences -Bulgaria) (Figure 3). LINKAM Scientific Instruments THMS600 provided with T95-Linksys PC interface and Linksys 32X system controller software, enabling PC control of temperature stability <0.1°C, gas flow (linear cooling speeds from 0.01 to 100°C/min), data acquisition and export, as well as multiple rise/down programming. The sample area is 22 mm diameter. A 100 Ohm platinum resistor sensor makes available more accurate and stable temperature signal that can be achieved with a thermocouple.

PixeLINK application software we use for digital image capture of the PCM slurry or PCM composite morphology in parallel to their thermal behaviour:  Temperature range of phase transition, Tonset, Tmax compared to DSC data.

 Crystallization/nucleation phenomena in situ.

 Stability of interfacial surfactant molecular layer in PCM emulsion slurry.

 Droplet size, shape, distribution.

 Temperature effect on the salt adhesion process in the graphite/Al matrix.

Scanning electron microscopy

Scanning electron microscopy is a useful tool to evaluate morphological characteristics of PCM and PCM composites at the microscale. It has been used by CERTES -France & Polymer Institute (France) and the Slovak Academy of Sciences (Slovakia) to evaluate PCM particles shape, the quality of the silver coating PCM and also the dispersion of fillers inside the polymeric matrix. SEM (Jeol 6301F) fitted with an X-ray energy dispersive spectrometer (EDS, Oxford, Link Isis) was used and SEM images like the one presented in Figure 4 were taken at 0° tilt angle, using a 10-15 kV accelerating voltage. The samples (PCM particles and PCM particles metalized with silver and composites) were covered previously with a thin layer of gold to guarantee their electrical conductivity in order to prevent drift and charging.

At Riga Technical University (Latvia), a simple method was used for preparing samples for determination of the shell thickness of PCM microcapsules. The slurry of microcapsules was dried in the open air. The dried specimen was broken and fracture was exposed to SEM observation (Figure 5). Riga Technical University also used SEM to investigate damage resistance of microcapsules in the cement-lime mixtures before and after thermal cycling. For example, Figure 6 shows that the microcapsules were not damaged by the mixing procedure of mortar and cycling for 300 times in the temperature range from +5 to 45ºC [START_REF] Bajare | Development of cement and lime based plaster with microencapsulated Phase Change Materials (PCM)[END_REF].

Atomic Force Microscopy

At the Institute of Physical Chemistry Romanian Academy (Romania), in the study of PEG composites, AFM is used to evaluate the surface morphology of the P15-E. AFM Dynamic Force Module (intermittent contact mode), using an EasyScan 2 apparatus (Nanosurf AG, Switzerland) with a high resolution scanner (10 x 10 μm with vertical range of 2 μm and z-axis resolution of 0.027 nm) and a X-Y linearity mean error less than 0.6%. The scan rate was about 1-2 Hz. The phase contrast working mode for imaging employed sharp Si tips (NCLR from Nanosensors™) with a radius of curvature less than 10 nm (typically 7 nm). The 1.5x1.5 μm phase image of the P15-E in the Figure 7 shows a lamellar morphology with an interlamellar spacing of 55 nm. The XRD calculation of the d-spacing for the 120 reflection shows is tantamount to 0.52784 nm.

At the University of Barcelona (Spain) AFM is used to evaluate the morphology and mechanical integrity of microencapsulated PCM with temperature [START_REF] Giro-Paloma | Stiffness characterization of microencapsulated phase change materials by atomic force microscopy[END_REF]. The aim of this study is to evaluate the mechanical performance of microencapsulated Micronal ® DS 5001 provided by BASF ® at different temperatures of use: 25 ºC, 45 ºC and 80 ºC, approximately. When microencapsulated PCM is used as aqueous slurry in active storage systems, changes are observed after several thermal cycles that are attributed to a partial degradation of the microcapsules by breakage. For this reason, it is interesting to measure the maximum force these particles can hold before break.

Observation of Micronal particles with SEM reveals that small particles (6 μm diameter) aggregate in bigger particles around 150 μm). Both types of particles were indented with the cantilever using the following experimental conditions:

 Aggregates (100M) 150 μm F = 3 μN  Single particles (10M) 6 μm F = 1.5 μN
In Figure 8 results of the several indentations are represented statistically as the frequency of each outcome for aggregates and individual Micronal particles at different temperatures.

Molecular spectroscopies

Raman spectroscopy is a very useful technique in studying chemical structure and morphology of heterogeneous materials as the PEG composites. At the Institute Physical Chemistry of the Romanian Academy (Romania) unpollarized Raman spectra of the P15-E composite at room temperature were obtained on a Jasco NRS-3300 micro-Raman spectrometer with an air cooled detector in a Backscattering geometry using a microscope objective of 100x. Two excitation laser lines, operating at 488 and 785 nm, were used with the power at the sample surface of 12.3 and 5.2 mW, respectively. Spectra were recorded in a single subtractive configuration (gratings 1800 gr/mm and 600 l/mm) over a wide frequency range of 100-4000 cm -1 for Ar + -ion laser and 100-1800 cm -1 for diode laser. The fluorescent background for each spectrum of the P15-E composite was subtracted by using linear segments. All corrected Raman spectra were fitted with Lorenzian profiles by means of Igor software. Unpolarized Raman spectrum of the P15-E composite and its spectral deconvolution is shown in Figure 9 [START_REF] Constantinescu | Latent heat nano composite building materials[END_REF].

Unpolarized solid state Raman spectra were recorded over 10-900 cm -1 and 200-1800 cm -1 using a RM-1000 RENISHAW Raman Microscope, equipped with a charge coupled device (CCD) detector, and a Leica microscope with a 50x objective, in backscattering geometry. The emission line at 514.532 nm of an Ar + laser was used to excite the samples. The diameter of the laser spot on the sample surface amounted to 2-3 m providing a spectral resolution better than 2 cm -1 . A THMS-600 cell (LINKAM) was used for temperature control of the samples from 25-600ºC with an accuracy of the 0.1ºC. Raman spectra collected at different temperatures were temperature-reduced to account the first order Bose-Einstein distribution factor from the experimental spectra. As is shown in Figure 10, the images and band at 2500-3500cm -1 indicate the influence of temperature at the phase change by cycling.

At the University of Barcelona (Spain) Fourier Transformed Infrared Spectroscopy (FT-IR)

with an ATR accessory was used to characterise wood-PCM composites prepared at Riga Technical University (Latvia) [START_REF] Vecstaudza | Characterization of wood/PCM composites for building applications[END_REF]. In this study, black alder wood (Alnus Glutinosa) was impregnated with paraffins RT21 and RT27 -PCMs manufactured from Rubitherm GmbH.

Samples were also coated with a polystyrene based coating to avoid leakage. Introduction of PCM in wooden matrix gives rise to absorption bands in the wavenumber region of 3000 -2800 cm -1 . These bands correspond to methyl and methylene group stretching, as RT21 and RT27 are saturated hydrocarbons with molecular formula C n H 2n+2 . Wood adds characteristic absorption bands of cellulose, hemicellulose and lignin to the wood/PCM composite spectra. Bands of polystyrene usually appear of composite layers closer to surface or ends of the sample.

Homogeneity of impregnated PCM was evaluated by comparing spectra of wood/PCM composite cut in several slices. Differences between slices of one sample are illustrated in Figure 11.

At Çukurova University (Turkey), PCM samples are subjected to melting/freezing thermal cycles carried out to determine their thermal and chemical stability. The number of thermal cycles varies depending on the application of PCM. FTIR analysis is used to test the chemical stability. The basic approach is to carry out FTIR analysis before and after thermal cycling and to compare the characteristic peaks obtained in the spectra. Figure 12 shows FT-IR spectra of the non-cycled and cycled (200, 400, 600 cycles) Capric acid-Myristic acid (CA-MA) binary mixture as PCM [START_REF] Cellat | Enhancing Thermal Properties of Concrete Mixtures with Phase Change Materials, 2nd International Conference on Sustainable Energy Storage in Buildings[END_REF]. As seen in the figure, all peaks fit into one another at the same frequency band. This means that repeated thermal cycles did not cause any degradation in the chemical structure of PCM the mixture.

Physico-chemical

Chemical stability related to fire testing

For testing of PEG-epoxy composite P15-E, a modified set up in comparison to the ASTM E84

Standard Test Method for Surface Burning Characteristics of Building Materials was used at the Institute of Physical Chemistry of the Romanian Academy (Romania), where a reduced size specimen with 150 mm long, 75 mm wide and 25 mm thick was placed between two plaster sheets of 125 mm thickness. To avoid direct contact of the fire with the sample the two plaster panels were 20 mm longer than the rectangular specimen, as is shown in Figure 14 left. Thus, indirect slow heating at the bottom side of the material cuboid was used. The temperature on both sides of the sample was monitored by means of two K thermocouples in time (Figure 14 right). The bottom side temperature of the specimen did not exceed 100ºC during the 20 minute span of the test in order to monitor the thermal degradation of the PCM component, i.e. PEG 1500, leaving unaffected the thermosetting component (the cured epoxy resin) prone to degradation at higher temperatures.

Degradation of PCM composites after the fire testing was evaluated using FTIR with coupled ATR Spectra of the FP15-E and RP15-E samples along with their Gaussian deconvolutions within 1200-975 cm -1 and 3700-2600 cm -1 ranges. The least varying peak at 1467 cm -1 (CH 2 deformation) was chosen for spectra normalization of the spectra. The IR spectra emphasize these changes produced inside the material.

Figure shows the Gaussian deconvoluted IR spectra of the RP15-E and FP15-E samples in the 1200-2400 cm -1 range, where the peaks of ester groups appear after the fire testing, that can be attributed to the thermo-oxidative degradation of the PEG. Chemical compatibility between the PCM, the shape stabilized PCM, PCM slurries and the materials containing them are usually carried out following protocols of corrosion tests (for example, at the University of Lleida and the University of Zaragoza, Spain). For inorganic PCM corrosion of metallic containers is given as a rate of corrosion. Results on this can be found in Cabeza et al. and Oró et al. [START_REF] Cabeza | Immersion corrosion tests on metal-salt hydrate pairs used for latent heat storage in the 32 to 36°C temperature range[END_REF][START_REF] Cabeza | Middle term immersion corrosion tests on metal salt hydrate pairs used for latent heat storage in the 32 to 36ºC temperature range[END_REF][START_REF] Cabeza | Immersion corrosion tests on metal-salt hydrate pairs used for latent heat storage in the 48 to 58ºC temperature range[END_REF][START_REF] Cabeza | Long term immersion test son metal-PCM pairs used for latent heat storage in the 24 to 29 ºC temperature range[END_REF][START_REF] Oró | Corrosion of metal and polymer container for use in PCM cold storage[END_REF]. Otherwise when plastic containers are considered, migration of PCM within the container and water sorption through the plastic wall is usually evaluated by gravimetric analysis after cycling the PCM. Lázaro et al. [START_REF] Lázaro | Experimental study on phase change materials and plastics compatibility[END_REF], Castellón et al. [START_REF] Castellón | Compatibility of plastic with pahse change materials (PCM)[END_REF], and Oró et al. [START_REF] Oró | Corrosion of metal and polymer container for use in PCM cold storage[END_REF] reported some studies of compatibility of plastics with PCM.

Compatibility between PCM and container materials

The methodology used includes that once the samples were weighted, they were immersed inside glass test tubes containing the PCM formulations. Then, the glass tubes were covered with plastic paraffin to avoid contact with environmental agents. The evolution of corrosion rate in time was evaluated. Visual evaluation was done seeking for bubbles, precipitates, surface changes and pitting process. pH measurements with litmus paper were also carried out. Before re-weighted, samples were cleaned and dried. While the metal samples were polished with abrasive paper and then dried, the polymer samples were cleaned with soft paper and dried with compressed air.

Results of compatibility of PCM with metals are given on a manufacturers recommendation (Table 1) and are shown in Table 2. Figure 16 shows some samples after the experimentation.

On the other hand, in order to evaluate the polymeric materials, mass change Δm(%) was calculated. The mass changes could be positive or negative, depending if the polymer absorbs or is diluted by the PCM, respectively.

Results show that polymers here analysed presented no significant mass gain (Table 3). Hence, those materials are suitable and recommended to be used as PCM container with the PCM studied in these papers. Polymers are materials commonly used as containers in other application due to their low price, low density and easy shaping.

Rheological behaviour

A TA Instruments AR-G2 device (a controlled stress rheometer) is used at the University of Zaragoza (Spain) to study the rheological behaviour of the materials [START_REF] Delgado | Determining the rheological behavior of octadecane as phase change material: First approach[END_REF]. Depending on the objective, different studies can be developed. For those materials that will be pumped in real applications, the viscosity vs. share rates and the viscosity vs. temperature curves are obtained.

For those materials that will remain at rest, viscosity vs. temperature curves are obtained with share rates close to zero. These last measurements can be complemented with measurements made in oscillatory mode. Figure 17 shows examples of these results obtained for octadecane at the University of Zaragoza.

Thermal stability investigation

A calorimeter was designed at the Warsaw University of Technology (Poland) for the performance of fast melting and solidification (cycling) of phase change materials and investigation of the stability of their thermophysical properties [START_REF] Domanski | Measurements of thermal properties of PCM materials[END_REF]. Since the volume of sample container equals to 200 cm 3 , results are representative for real elements with PCM which are used in thermal storage units. Temperature range is limited to 20-150C, mainly because the cooling system used in the set-up.

The scheme of the main part of the calorimeter is shown in Figure 18, and it is used to test the thermal stability of PCM.

Cycling tests consists of two consecutive phases, which are repeated from several to several hundred times, depending of the thermal behaviour of the material. In order to check whether material is stable or not, temporal temperature variations during melting and solidification are recorded and analysed. In Figure 19 temperature vs. time for selected cycles of the tests performed on the hydrate Na 2 S 2 O 3 5H 2 O are presented. Since the temperature characteristics change substantially in subsequent cycles, this material is not stable. In 28 th cycle melting is not observed, that means that the hydrate decomposed into the salt and water. If material is stable all recorded characteristics look very similar to that for the first cycle.

Physical stability testing of slurries

Microencapsulated paraffin based phase change slurry (PCS) from BASF has been characterized at the Fraunhofer Institute for Solar Energy Systems (Germany). The initial slurry containing 40 %wt. of capsules was diluted to a capsule fraction of 20 %wt. by adding water.

As the density between microencapsulates and water is different, a thickener was added to slow down the separation speed. The enthalpy measured by Micro DSC III from Setaram reveals a melting range between 16 and 18.5 ºC, a melting enthalpy of about 40 kJ/K and a hysteresis for the solidification curve of 2 K.

After characterization the PCS was tested with a storage test facility. It contains plate type heat exchangers, centrifugal pumps, a 500 L storage tank and sensors to measure all relevant volume flows and temperatures (Figure 20). The heat which is transferred over the heat exchanger from water into the PCS is stored in the tank. Therefore, the volume flow rate through and the temperature difference over the heat exchanger on the thermostat side (red in the schematic drawing -Figure 20) was measured. The heat losses of the storage have been determined with previous experiments using water instead of PCS. For the measurements the storage is operated in different temperature intervals. In the first experiment the intervals have been kept always at 2 K and the temperature level was increased in 2 K steps starting from 12 °C (Figure 21 left). In a second experiment the temperature intervals have been increased from 2 K to 6 K and the stored heat was compared with that reached by using water within the same temperature intervals (Figure 21 right).

Mechanical properties

Mechanical characterization at nanometric scale

The grid indentation technique has been employed to isolate the mechanical properties of a composite material or graphite foam and salt at the University of Barcelona (Spain) in collaboration with TREFFLE laboratory (France). This technique is a powerful tool for examining the properties of constituent phases independently of each other in composite material microstructures without the need to observe the residual imprints. Moreover, the mechanical properties (hardness and Young's modulus) for individual phases can be isolated, as well as their surface volume fractions at nanometric scale. This technique and the statistical analysis by the deconvolution process allow representing the mechanical properties for each phase as a normal density function [START_REF] Canseco | Mechanical characterization at nanometric scale for heterogeneous graphite-salt phase change materials with a statistical approach[END_REF]. It can also be employed to other composite materials containing mixture of constituent phases distributed randomly. In Figure 22 is represented the Young's modulus distribution simulated using the CDF fitting parameters. To complete the mechanical characterization, a macro structure mechanical characterization like micro compression tests must be done.

Loss of mechanical properties of polymers in contact with paraffin. Evaluation using nanoindentation

Nanoindentation technique is adequate to evaluate mechanical changes in the surface and in the internal properties of polymers when exposed to paraffin for a long time. Moreover, the variation of the mechanical parameter with depth of indentation gives also valuable information that may be suitable for design purposes. In this study, carried out at the University of Barcelona (Spain) [START_REF] Giró-Paloma | Characterisation of thermoplastics with nanoindentation[END_REF], polypropylene PP unfilled and filled with Mg(OH) 2 (40 phr) were submerged in melted paraffin Rubitherm ® : RT-25 at different temperatures for 32 days.

Stiffness and hardness of polypropylene decreases due to contact with paraffin, and the effect is more significant with increasing temperature. Table 4 summarises some of the results.

Conclusions

This paper presents different methods to characterize structural and morphological properties of PCM [START_REF] Delgado | Determining the rheological behavior of octadecane as phase change material: First approach[END_REF].

rmal propert ties of

Stability 3 . 1 .

 31 ChemiluminiscenceChemiluminiscence results from the conversion of chemical energy into light which takes place in many elementary reactions of more complex reaction mechanisms. During oxidation the molecules of a reaction product which are formed in an excited state are either carbonyl compounds or singlet oxygen. These products are formed during recombination of peroxy radicals.The equipment unifies knowledge acquired in the study of chemiluminescence of organic substrates since 1960 by groups of researchers in Australia, Russia, United Kingdom, USA and the former Czechoslovakia (Polymer Institute of Slovak Academy of Sciences). It is suitable for the detection and measuring of thermal oxidation of polymers or thermal decomposition of preoxidized polymers containing hydroperoxides from ambient temperature to 240 °C. The standard atmosphere is ensured by flow of gas at maximum rate of 150 ml/min. Extreme sensitivity, which is achieved by almost zero background noise, enables to measure the oxidation of polymers already at 40 °C and the tendency towards oxidation of organic substrates even in the so-called "induction period" may be well estimated.The final form for the chemiluminescence intensity is expressed as the sum of two processes composing the resulting light emission. The degradation of the PEG composites developed at the Institute Physical Chemistry of the Romanian Academy (Romania)[START_REF] Constantinescu | Latent heat nano composite building materials[END_REF], was studied and the two processes are related to the rate of initiation of polyethylene glycol (PEG) degradation via pre-exponential factors and activation energies. All the kinetic parameters, i.e. pre-exponential factor, activation energy, and average kinetic rate, are calculated for 150 °C. Modelling of the kinetics is shown in Figure13.

  Figur

  Figur

8 Figure 22 .

 822 Figur coolin coolin

Acknowledgements

The work is partially funded by the European Union (COST Action COST TU0802) and the Spanish government (ENE2011-28269-C03-01, ENE2011-28269-C03-02 and ENE2011-22722). The authors would like to thank the Catalan Government for the quality accreditation given to their research group GREA (2014 SGR 123) and their research group DIOPMA (2014 SGR 1543). Aran Solé would like to thank the Departament d'Universitats, Recerca i Societat de la Informació de la Generalitat de Catalunya for her research fellowship.