Quy-Dong To
email: quy-dong.to@univ-eiffel.fr

Guy Bonnet

Fourier Transform approach to homogenization of frequency dependent heat transfer in porous media

Keywords: Fourier Transform method, Strong contrast homogenization, Heat transfer, Dynamic thermal permeability, Porous media

The Fourier Transform method is used to solve frequency dependent heat transfer problems. The periodic Lippmann-Schwinger (LS) integral equation in Fourier space with source term is first formulated using discrete Green operators and modified wavevectors which can then be solved by iteration schemes. The objective is to solve the local problem involving strong contrast heterogeneous conductive material, with application to gas-filled porous media with both perfect and imperfect Kapitza boundary conditions at the bi-material interface. The effective parameters such as the dynamic conductivity and the thermal permeability in the acoustics of porous media are also derived from the cell solution. Numerical examples show that the schemes converge fast and yield accurate results when compared with analytical solution for benchmark problems. The formulation of the method is built using static or dynamic Green operators and can be applied to pixelized microstructure issued from tomography images.

Introduction

Determining the overall behavior of heterogeneous media is the main focus of homogenization theory [START_REF] Milton | The theory of composites[END_REF][START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF][START_REF] Nemat-Nasser | Micromechanics: overall properties of heterogeneous materials[END_REF]. Based on the local boundary value problem on a Representative Volume Element (RVE), average quantities can be obtained to establish the macroscoscopic relations for the equivalent homogenized material. By construction, the distribution of the heterogeneous materials in the RVE can be either random or periodic.

The homogenization of periodic materials can be addressed by double scale asymptotic procedures [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF]. By developing the local problem with respect to the scale separation ratio ε, a large number of macroscopic models can be derived, whose form can be different from the local constitutive behavior. Considering the heat transfer problem under harmonic excitation at a given frequency ω, [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF] showed that for large contrasts of thermal conductivity, the constitutive equation is frequency dependent, resulting in a memory kernel for the transient behavior.

It is worthwhile noticing that the local cell problem is of the same type as the one related to "dynamic thermal permeability" that appears in the computation of the dynamic compressibility of a gas within a porous medium [START_REF] Lafarge | Dynamic compressibility of air in porous structures at audible frequencies[END_REF].

The main objective of the present work is to address the frequency dependent heat transfer problem in a periodic cell. With the increasing capabilities of computing power, solving problems on realistic microstructures becomes more accessible. Numerical methods are developed to handle both tomographic or computer generated images. Compared with other methods such as Finite Element Method or Random Walk, methods based on Fourier Transform benefit from serious advantages. The methods rely directly on pixelization of the simulation domain and the robustness of the Fast Fourier Transform algorithm when dealing with periodicity. Additionally, the resolution based on iterative schemes is fast and yields accurate overall properties.

Although the Lippmann-Schwinger (LS) integral equations had been derived previously [START_REF] Brown | Solid mixture permittivities[END_REF][START_REF] Kröner | Bounds for effective elastic moduli of disordered materials[END_REF] for micromechanics problems, [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF]; [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF] proposed to use iterative schemes combined with Fast Fourier Transform to solve numerically these equations. A broad range of engineering problems such as fluid flows [START_REF] Nguyen | A Fourier based numerical method for computing the dynamic permeability of periodic porous media[END_REF][START_REF] Bignonnet | Efficient FFT-based upscaling of the permeability of porous media discretized on uniform grids with estimation of RVE size[END_REF], non linear behavior [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF], finite strain [START_REF] De Geus | Finite strain FFT-based non-linear solvers made simple[END_REF], Kapitza interface problem in composites [START_REF] Monchiet | FFT based iterative schemes for composites conductors with non-overlapping fibers and Kapitza interface resistance[END_REF] etc.... have been treated by the method since then. Efforts have been made to accelerate the convergence at high contrast ratio such as the polarization schemes [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF][START_REF] Monchiet | A polarization-based fast numerical method for computing the effective conductivity of composites[END_REF], variational schemes [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF], Newton-Krylov solvers [START_REF] Brisard | Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFTbased numerical method for the homogenization of composites[END_REF][START_REF] Vondřejc | An FFT-based Galerkin method for homogenization of periodic media[END_REF][START_REF] Zeman | Accelerating a FFTbased solver for numerical homogenization of periodic media by conjugate gradients[END_REF][START_REF] Schneider | On the Barzilai-Borwein basic scheme in FFT-based computational homogenization[END_REF].

For porous materials where the contrast ratio becomes infinite, the convergence becomes the main issue. It is shown that most of iterative schemes fail to satisfy the divergence free and rotational free criteria [START_REF] Moulinec | Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials[END_REF]. To overcome those limitations, [START_REF] To | FFT based numerical homogenization method for porous conductive materials[END_REF]; [START_REF] To | Fourier transform approach to non periodic boundary value problems in porous conductive media[END_REF] recently proposed a class of iteration scheme based on boundary integral equation. Another interesting technique is to use discrete Green tensors and modified wavevectors as derived from finite difference approximation [START_REF] Willot | Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF]. In a comparison between different iterative schemes, [START_REF] Moulinec | Convergence of iterative methods based on Neumann series for composite materials: theory and practice[END_REF] the authors stated "continuous Green's operator performs better for moderate contrast whereas discrete operators should be used for large contrast". Since the latter is suitable for pixelized microstructure, it will be adopted to solve the heat transfer problem.

The paper is organized as follows. First, the important results of asymptotic homogenization theory in the case of large contrasts between the properties of constituents is recalled. Next, the LS integral equation for the heat transfer problem with a source term using modified Green tensors is presented. Combined with the immersed interface technique [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF], the iterative schemes are then derived for a large class of boundary value problems: prescribed flux, prescribed temperature, temperature jump, temperature dependent source term. These schemes are then used to solve the frequency dependent heat transfer problem and results are compared with analytical solutions from the literature.

Heat transfer problems from asymptotic homogenization theory

Strong contrast composites and effective frequency dependent heat capacity

According to the asymptotic homogenization theory, the macroscopic law is constructed based on the length scale ratio

ε = l micro L macro ≪ 1 (1)
where l micro and L macro are the characteristic length at the micro and macro scales of the heterogeneous material. In addition, all physical quantities are treated as functions of two space variables, the macroscopic space variable X and the microscopic space variable x, with x = X/ε representing local fluctuations within the microstructure.

Let us consider a porous material whose solid skeleton is characterized by conductivity κ s and thermal capacity c s filled by a fluid characterized by κ g and c g as in Fig. 1. When κ g /κ s = O(1) and ρ s c s /(ρ g c g ) = O(1) the asymptotic expansion method [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF] shows that the porous material behaves as a material with the effective conductivity given by the usual homogenization method, i.e the static case. In this case, the effective thermal capacity is equal to the volume average of local thermal capacities.

When the thermal conductivity of the gas is very small compared to the one of the solid, this result is no longer valid. It is shown that if κ g /κ s = O(ε 2 ) and ρ s c s /(ρ g c g ) = O(1), the porous medium is characterized by a frequency dependent thermal capacity that is obtained from a dynamic localization problem which writes:

κ g ∆ x θ -iωρ g c g θ -iωρ g c g T = 0, (2)
and θ = 0 at the solid/gas boundary (3)

where T (X) is the amplitude of the macroscopic temperature, that depends only on the macroscopic space variable X. The solution of this equation leads to a local value of θ linearly related to T :

θ(X, x, ω) = -χ(x, ω)T (X) (4)
The macroscopic balance of energy can be written:

div X (κ ef f ∇ X T ) = (ρc) ef f iωT (5)
where κ ef f is the effective conductivity and (ρc) ef f is the effective dynamic thermal capacity given by (ρc

) ef f = ⟨ρc⟩ - ρ g c g |V | Vg χ(x, ω)dx. (6)
Here, V is the volume of the periodic cell, V g is the volume of fluid within the cell and ⟨⟩ denotes the volume average over the cell.

Heat transfer due to acoustic wave in fluid filled porous media

In this case, the dynamic behavior is governed at the local scale by mass conservation, Navier-Stokes equation, heat equation and the equation of state for the gas. The scaling of these equations with suitable renormalization of the physical coefficients has been described in [START_REF] Auriault | Homogenization of coupled phenomena in heterogenous media[END_REF] and [START_REF] Boutin | Acoustic absorption of porous surfacing with dual porosity[END_REF]. The asymptotic method leads to the following results:

• At highest order in ε, the macroscopic pressure P depends only on the macroscopic space variable.

• At the following order, the local temperature within the gas is the solution of

κ g ∆ x θ -iωρ g c g θ + iωP = 0, (7)
with θ = 0 on the boundary of the fluid domain. It is noteworthy that equation ( 7) has the same form as (2) where the driving force iωP plays the same role as -iωρ g c g T

The solution is proportional to iωP and it is natural to introduce a variable θ given by θ = θκg iωP that is solution of:

∆θ - θ δ 2 + 1 = 0, ( 8
)
where δ is the complex length

δ = κ g iωρ g c g . ( 9
)
The volume average of θ over the fluid volume is given by:

Θ = 1 |V g | Vg θdx (10)
with Θ being a macroscopic physical parameter called the "thermal permeability" [START_REF] Lafarge | Dynamic compressibility of air in porous structures at audible frequencies[END_REF]. For a given microstructure, Θ only depends on the complex length δ, which is a function of frequency by (9), or Θ = Θ(δ(ω)). This parameter plays an important role in the determination of the effective dynamic compressibility.

From equations ( 2) and ( 7), it can be seen that these equations are the same when replacing P by -ρ g c g T . It follows that χ in (4) is related to θ * by

χ = θ δ 2 . ( 11
)
Finally, the effective dynamic capacity is obtained from the thermal permeability by:

(ρc) ef f = ⟨ρc⟩ - 1 δ 2 ρ g c g |V | |V g | Θ = ⟨ρc⟩ - 1 δ 2 ν g ρ g c g Θ (12)
where ν g is the volume concentration of gas.

In the following, the problem of dynamic thermal permeability will be solved, knowing that the dynamic thermal capacity can be obtained from this relation.

Temperature jump boundary conditions

In addition to the perfect boundary conditions (3), we also study the case where there is a temperature jump, also called temperature slip, in the form [START_REF] Kapitza | The study of heat transfert in helium II[END_REF]:

θ(x) = l∇θ(x)n(x) = κ g κ int ∇θ(x)n(x) at solid/gas boundary (13)
with n being the normal vector, l the slip length, κ int the interface coefficient. The phenomenon is frequently encountered at gas flows characterized by high Knudsen number, i.e when the mean free path is comparable to the characteristic size of the pore [START_REF] Umnova | Effect of boundary slip on the acoustical properties of microfibrous materials[END_REF][START_REF] Chastanet | Acoustics with wall-slip flow of gas-saturated porous media[END_REF]. Microscopically, it can be quantified by the collision between the gas atoms and the solid wall (Liao et al., 2018a,b). For perfect interface, l = 0 and κ int = ∞, the perfect boundary conditions (3) is recovered.

It is noteworthy that the slip condition can be recovered as the limit case of a layer of large conductivity and small thickness at the interface between the phases, when the thickness tends to zero. This does not affect the distribution of heat capacity and it will be considered that the result of homogenization theory [START_REF] Auriault | Homogenization of coupled phenomena in heterogenous media[END_REF] is also valid in this case, the difference being only to account for the temperature jump in the localization problem.

Fourier Transform resolution method using the static Green tensor

For simple geometries (plane slits, cylindrical interfaces), the boundary value problem (8-9) can be solved in a closed form, but in the general case a suitable numerical solution becomes necessary. In this section, a numerical solution based on a FFT formulation is presented.

Periodic heat equation with source term and associated discrete operators

Let us consider a heterogeneous material where the local conductivity κ(x) is a periodic function of the coordinate x with period V of dimensions V = L 1 × L 2 × L 3 in the Cartesian coordinate system. The governing equations of the problem are the following

e = E -∇θ, div(j) + s = 0, j = κe (14)
where θ, e, j and s are the V periodic temperature, (minus) temperature gradient, heat flux and source term with volume averages Θ, E, J and 0, i.e

⟨θ(x)⟩ = Θ, ⟨e(x)⟩ = E, ⟨j(x)⟩ = J , ⟨s(x)⟩ = 0. ( 15
)
The last condition ⟨s(x)⟩ = 0 is necessary for the compatibility between the periodic flux j and the energy conservation. In general, the source term s can be temperature dependent. Due to the periodicity, the problem can be solved using Fourier series. For a V periodic function f (x) of real variable x(x 1 , x 2 , x 3 ), it can be written as

f (x) = ξ f (ξ)e iξx , i = √ -1 (16) where f (ξ) is the Fourier transform of f (x), notation f (ξ) = F[f (x)]
, and function of wave vector ξ(ξ 1 , ξ 2 , ξ 3 )

ξ k = 2πn i /L i , n i = 0, ±1, ±2, ..., ±∞, i = 1, 2, 3. (17)
Numerically, we are limited to finite resolution by bounding n i with resolution parameter N , i.e -N ≤ n i ≤ N -1, and solve the problem on collocated grid. It is noted that the resolution parameter N corresponds to resolution 2N in real space. For example, if N = 128 is adopted for 2D problems, in real space the resolution would be 256 × 256.

The usual way to produce Green tensor for the FFT method is to use the classical relation between spatial derivative and the product by the wave vector. An alternative is to use a modified Green tensor built as follows using modified wave vectors. As shown in [START_REF] Willot | Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF], the use of modified Green tensors avoids the convergence issues in infinite contrast situations and solution fluctuations due to Fourier series. The spatial derivatives at the grid points x(n 1 h 1 , n 2 h 2 , n 3 h 3 ) with h i = L i /2N are estimated by forward and backward finite difference. In physical and Fourier space, we have Forward:

∂f ∂x i (x) = 1 h i [f (x + h i e i ) -f (x)] , ∂f ∂x i (ξ) = ik i f (ξ), Backward: ∂f ∂x i (x) = 1 h i [f (x) -f (x -h i e i )] , ∂f ∂x i (ξ) = ik * i f (ξ), i = 1, 2, 3. ( 18
)
with the modified complex wavevector k i and its complex conjugate k * i being the quantities

k i = 1 h i (e iξ i h i -1), k * i = 1 h i (e -iξ i h i -1). ( 19
)
Next, following [START_REF] Willot | Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF], the discrete Green operators P , Q, H and R are constructed from the modified wavevector k(k 1 , k 2 , k 3 )

P (ξ) = k ⊗ k * |k| 2 , Q(ξ) = I -P (ξ), H(ξ) = I -2P (ξ) R(ξ) = ik |k| 2 , ∀ξ ̸ = 0 (20)
and

P (0) = Q(0) = H(0) = R(0) = 0 (21)
As a result, the integral equations to solve the heat transfer problem with source term can be obtained [START_REF] To | Fourier transform approach to non periodic boundary value problems in porous conductive media[END_REF]. Depending on the choice of unknowns, we have, say, for gradient e e(x)

= E + 1 κ 0 R * s -P * κ -κ 0 κ 0 e (22)
for the flux j

j(x) = J + R * s + Q * κ -κ 0 κ j (23)
and for the polarization τ

τ (x) = T + 2R * s + H * κ -κ 0 κ + κ 0 τ , (24)
Here the polarization τ with average ⟨τ (x)⟩ = T , is given by

τ = (κ + κ 0 )e = (κ + κ 0 )j/κ (25)
and * represents convolution operation in the real space. In Fourier space, this operation is equivalent to simple product between the transformed quan-tities.

The LS equations ( 22), ( 23) and ( 24) can be used to solve the corresponding quantities e, j, τ by iteration schemes. In addition, the temperature distribution can be recovered from the equivalent relation of ( 14)

θ(x) = Θ + θ(x), θ(x) = -R * * e. ( 26
)
and θ(x) is the fluctuation θ(x) around the mean value Θ.

Iteration schemes

Let us start with the LS equations ( 22), ( 23) and ( 24). The common points of those equations are that they have the form

u(x) = U + A * s + B * (ϵu) (27)
where u (which can be either e or j or τ ) represents the unknown to be obtained, U the volume average of u(x), ϵ(x) the known functions of the local conductivity, A, B are discrete Green operator with an explicit expression in Fourier space (see Table 1).

In section 2, it has been shown that we need to solve equations ( 8) with a temperature dependent source term. We consider the case where the source term s is a linear function of the local values of e or j or θ, i.e

s(x) = α e (x)e(x) + α j (x)j(x) + α θ (x)θ(x) + α 0 (x) (28)
where α e (x), α j (x), α θ (x) and α 0 (x) are known fields, thus independent from the solution e or j or θ. Since the source term must have zero mean ⟨s(x)⟩ = 0 and θ(x) = Θ + θ(x) from ( 26), the following relation must be verified The problem ( 14) can be solved using any one of the four governing equations ( 22), ( 23) or ( 24) and iterate in order to produce the field θ that is consistent with the dynamic equation. Algorithm 1 shows the case of polarization scheme, i.e u = τ that is used in the paper. The convergence quality of polarization method has been demonstrated in previous publications (see e.g [START_REF] Willot | Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF][START_REF] Moulinec | Convergence of iterative methods based on Neumann series for composite materials: theory and practice[END_REF][START_REF] Monchiet | A polarization-based fast numerical method for computing the effective conductivity of composites[END_REF].

-⟨α θ (x)⟩Θ = ⟨α e (x)e(x) + α j (x)j(x) + α θ (x) θ(x) + α 0 (x)⟩. (29) Name u(x) U A(ξ) B(ξ) ϵ(x) e scheme e(x) E R(ξ)/κ 0 -P (ξ) κ(x)-κ 0 κ 0 j scheme j(x) J R(ξ) Q(ξ) κ(x)-κ 0 κ(x) τ scheme τ (x) T 2R(ξ) H(ξ) κ(x)-κ 0 κ(x)+κ 0
Algorithm 1 Schemes for periodic heat transfer problem with source term Initialization τ (0) = T repeat 1.

Compute e (n-1) = τ (n-1) κ(x)+κ 0 , j (n-1) = κ(x)τ (n-1) κ(x)+κ 0 , θ (n-1) = -F -1 {R * F[e (n-1) ]} 2. Compute Θ (n-1) from e (n-1) , j (n-1) and θ (n-1) by ( 29) and θ (n-1) = θ (n-1) + Θ (n-1) 3. Compute s (n-1) from θ (n-1) , j (n-1) , e (n-1) by ( 28

) 4. Update τ (n) = T + F -1 {2RF[s (n-1) ] + HF[ κ(x)-κ 0 κ(x)+κ 0 τ (n-1) ]}. Equiva- lently, τ (n) = τ (n-1) -2F -1 {P F[j (n-1) ] -RF[s (n-1) ] + κ 0 QF[e (n-1) ]} n ← n + 1 until convergence criteria is met
Some remarks on the iterative schemes dealing with periodic problems can be made as follows.

i) In the classical homogenization problem where s = 0 and with a finite contrast, the local conductivity being bounded by κ min < κ(x) < κ max , the best choice of the reference conductivity in this case is κ 0 = √ κ min κ max for τ scheme. The polarization τ scheme has shown to have a faster convergence rate than e and j schemes. Furthermore, combined with the modified Green operators, it does not encounter infinite contrast issues.

ii) In the case where the source term s is prescribed s = α 0 (x) with the condition ⟨α 0 (x)⟩ = 0, step 3 is unnecessary. The temperature field θ is defined up to a constant Θ and the constant Θ (n-1) in step 2 can be set unchanged and equal to Θ during the whole process. In the opposite case where s depends on the local temperature θ, the constant Θ (n-1) in step 2 is then chosen to guarantee ⟨s (n-1) ⟩ = 0 by ( 29). See also applications to Neumann and Dirichlet boundary conditions in iii) and iv). The method can be extended to non periodic problems where the source term can be used to impose heat flux and temperature, for example the immersed interface technique [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF][START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF][START_REF] Wiegmann | Fast Poisson, fast Helmholtz and fast linear elastostatic solvers on rectangular parallelepipeds[END_REF]. The details of FT implementation using form factor and continuous Green tensors can be found in a recent work [START_REF] To | Fourier transform approach to non periodic boundary value problems in porous conductive media[END_REF]. Another possibility is to use polarization (or eigenstress/strain for elasticity) as done by [START_REF] Gelebart | A modified FFT-based solver for the mechanical simulation of heterogeneous materials with Dirichlet boundary conditions[END_REF] or [START_REF] Chen | Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron x-ray micro-computed tomography[END_REF]. For the present problem, the modified Green tensors [START_REF] Willot | Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF] as shown previously and the immersed interface method to deal with perfect or imperfect boundary conditions on pixelized or curved boundary are combined.

iii) Considering the Neumann boundary conditions for normal flux j n (x)

j n (x) = j p n (x), x ∈ Γ j ( 30
)
on boundary Γ j which is adjacent to the void phase. Here j p n (x) is the prescribed normal flux. The void can be considered as a material with zero conductivity. If the Γ j is approximated to follow the pixel boundary, crossing the cells located at x 1 , x 2 , ..., x M (see Fig. 2), the normal flux at x i is related to the source term at x i by

s(x i )h = j p n (x i ) ( 31
)
and h is the cell dimensions. In the general case where Γ j is a known curve boundary curvature significantly different from the previous pixel based boundary, the source term can be evaluated by a more general formula

h 3 s(x i ) = Γ i j j p n (x)ds ≃ j p n (x i )|Γ i j | (32)
where Γ i j is the portion of Γ j in the cell i.

iv) The Dirichlet type boundary conditions for temperature

θ(x) = θ p (x), x ∈ Γ t (33)
are applied to a boundary Γ t and the source term is known outside Γ t . In this case, θ p (x) is the prescribed temperature. The boundary Γ t is not necessary adjacent to the void phase as Γ j of the previous case. Again, the source term s can be used to control the temperature in the cell crossed by Γ t , for example θ(x i ) = θ p i , i = 1, 2, .., M . Such constraint can be imposed by the penalty technique, leading to

s (n) (x i ) = κ pen (θ (n-1) (x i ) -θ p i ) ( 34
)
with large value of κ pen . The other technique is to use the augmented Lagrange technique, coupled with iteration

s (n) (x i ) = s (n-1) (x i ) + κ L (θ (n-1) (x i ) -θ p i ) ( 35
)
where κ L is the Lagrange coefficient, usually taken κ L = κ 0 . The integration constant Θ (n-1) in step 2 can be chosen to ensure the vanishing average source

x i ∈Γt s (n) (x i ) + x i / ∈Γt s(x i ) = 0 (36)
at each time step n by (29). Specifically, we have for penalty technique

x i ∈Γt κ pen ( θ (n-1) (x i ) + Θ (n-1) -θ p i ) + x i / ∈Γt s(x i ) = 0 (37) or Θ (n-1) = 1 M    x i ∈Γt κ pen (θ p i -θ (n-1) (x i )) - x i / ∈Γt s(x i )    (38)
and for augmented Lagrange technique

x i ∈Γt s (n-1) (x i ) + κ L (Θ (n-1) + θ (n-1) (x i ) -θ p i ) + x i / ∈Γt s(x i ) = 0 (39) or Θ (n-1) = 1 M    κ L (θ p i -θ (n-1) (x i )) - x i / ∈Γt s(x i ) - x i ∈Γt s (n-1) (x i )    . (40)
v) The temperature jump condition on the boundary Γ int adjacent to the void. The normal flux on Γ int is proportional to the temperature difference and κ int is the interface coefficient

j n (x) = κ int (θ(x) -θ p (x)), x ∈ Γ int . (41)
In this case, the heat exchange at x 1 , x 2 , ..., x M where the material is in contact with the temperature θ p i , i = 1, 2, .., M is given by

hs (n) (x i ) = κ int (θ (n-1) (x i ) -θ p i ) (42)
for a pixel based boundary and by

h 3 s (n) (x i ) = κ int (θ (n-1) (x i ) -θ p i )|Γ i int | (43)
for a curved boundary.

vi) The standard convergence criteria of the algorithm is based on the integrability of e

ε e = ∥ik(ξ) × e (n) (ξ)∥ F < ε e (44)
and the conservation of energy associated to j

ε j = ∥ik * (ξ) j (n) (ξ) + s (n) (ξ)∥ F < ε j (45)
which results in the combined criteria

ε je = max{ε e , ε j } < ε je = min{ε e , ε j } (46)
where ε e , ε j , ε je are error norms and ε e , ε j , ε je ≪ 1 are the smallness parameter controlling the accuracy of the solution. In the case where there are temperature constraints at x 1 , x 2 , ..., x M , an additional error tolerance ε θ on temperature is adopted

ε θ = max i=1..M |θ (n) (x i ) -θ p i | < ε θ . (47)
The global error and global tolerance are denoted as

ε = max{ε je , ε θ } < ε. (48)
The frequency dependent heat transfer problem is the particular case of the general problem considered in the previous section.Two phases are present: the gas (conductivity κ g ) and the solid (conductivity κ s ). The temperature to solve is called the excess temperature θ which is considered as being null (θ(x) = 0) in the solid phase. Due to the similar nature of the two equations ( 7) and ( 2), the acoustic wave problem (2) in gas filled porous material will be dealt as an example.

In this case, the average flux, gradient and polarization do not exist E = J = T = 0 and the internal field is induced by the source term s which is known in the gas phase s = -iωP + iωρ g c g θ in the gas phase (49) and unknown in the solid phase.

The two types of boundary conditions can be treated by Algorithm 1. For ideal boundary conditions, the augmented lagrangian with temperature constraints is applied to either the boundary cell or the whole solid phase. In this case, the value of the conductivity κ s does not affect the solution of our problem but its choice can have impacts on the convergence behavior. The numerical experiments show that setting κ s = 100 -1000κ g and global tolerance ε = 0.001 -0.0001 will yield the optimal value between computation cost and accuracy. Regarding imperfect boundary conditions, we set κ s = 0 so that the normal flux is equal to the source term as before. Additionally, since the resulting flux is proportional to the measure of the boundary and thus affects the source term s, results can be sensitive to the choice of pixel based or curved boundary.

Generalized scheme based on frequency dependent Green operators

Up to now, the Green operators P , Q, H and R are obtained for a homogeneous media of conductivity κ 0 and the heterogeneity effect of conductivity is accounted for by the polarization. For systems involving variation of internal heat energy, the source term s can be further decomposed into a heat generation s and a source term coming from the thermal capacity as

s = s + iωρcθ ( 50
)
where ρc is the heat capacity. Then, the Green operators can be generalized to account for the heat capacity of the reference media ρ 0 c 0 . It is noted that a similar treatment is used in the context of homogenization of Bloch wave in elastodynamics (see e.g Srivastava and Nemat-Nasser, 2012, and the references therein) For the present heat transfer problem, let us pose

s ′ = s + iω(ρc -ρ 0 c 0 )θ = s -iωρ 0 c 0 θ (51)
This shows that the use of the reference capacity introduces a new polarization term iω(ρc -ρ 0 c 0 )θ. Inserting this relation into ( 22) in Fourier space and noting that θ(ξ) = -R * (ξ) e(ξ) from ( 26) and R(ξ)R * (ξ) = 1/|k| 2 via (20), we obtain

1 + iωρ 0 c 0 κ 0 |k| 2 e(ξ) = 1 κ 0 R(ξ) s ′ (ξ) - 1 κ 0 P (ξ) (κ -κ 0 )e(ξ), ∀ξ ̸ = 0. ( 52
)
Passing the term 1 + iωρ 0 c 0 κ 0 |k| 2 to the right side of the equation and returning to the real space yields

e(x) = E + 1 κ 0 R ′ * s ′ - 1 κ 0 P ′ * (κ -κ 0 )e, ( 53
)
with the new Green tensors

R ′ , P ′ R ′ (ξ) = ik |k| 2 + iωρ 0 c 0 κ 0 , P ′ (ξ) = k ⊗ k * |k| 2 + iωρ 0 c 0 κ 0 , ∀ξ ̸ = 0. ( 54
)
Similar expressions can be derived for flux j and polarization τ j

(x) = J + R ′ * s ′ + Q ′ * κ -κ 0 κ j, (55)
and

τ (x) = T + 2R ′ * s ′ + H ′ * κ -κ 0 κ + κ 0 τ , (56)
with Q ′ , H ′ defined as

Q ′ (ξ) = I -P ′ (ξ), H ′ (ξ) = Q ′ (ξ) -P ′ (ξ), ∀ξ ̸ = 0 (57) Like R, P , Q, H, the tensors R ′ , P ′ , Q ′ , H ′ vanish at ξ = 0, or P ′ (0) = Q ′ (0) = H ′ (0) = 0. It is noted that if the composite ratio ωρ 0 c 0 /κ 0 → 0, we recover the original expressions R ′ = R, P ′ = P , Q ′ = Q, H ′ = H.
We find that LS equations ( 53), ( 55) and ( 56) based on the frequency and heat capacity dependent Green operators P ′ , Q ′ and H ′ have the same structure as ( 22), ( 23) and ( 24). As a result, the iteration method can be based on Algorithm 1 with suitable modification. The modified source term s ′ becomes

s ′ = -iωP + iω(ρc -ρ 0 c 0 )θ in the gas phase (58)
and is unknown in the solid phase. To enforce the given temperature at the solid phase, the penalty and the augmented Lagrange techniques are applied to s like in the original scheme to ensure (36) and s ′ is calculated by the formula

s ′(n) = s (n) -iωρ 0 c 0 θ (n-1) . ( 59
)
As a result when ρ 0 c 0 = 0, the scheme associated to the frequency dependent tensor is exactly the same as in the original scheme.

Numerical applications

In the previous sections, two kinds of Green tensors, the static and dynamic ones were presented. Since the first one leads to converging results, the results are presented in the first two subsections by using the static Green tensor . Next, the dynamic Green tensor is introduced in the last subsection in order to evaluate the capacity of the related numerical scheme, compared to the one related to static Green tensor.

Heat transfer with an ideal boundary

The two 2D problems where analytical solution is available (see Appendix) will be considered first. The unit cell has dimensions L 1 = L 2 = 1 and the resolution parameter is fixed at N = 128, or resolution 256 × 256 in real space. Other parameters of the simulation are: κ g = 1, κ s = 100, κ 0 = √ κ g κ s = 10 and global tolerance ε = 0.001. The ideal boundary condition is applied at the solid gas interface and the augmented Lagrange technique is used to control the temperature of the solid phase. The iteration scheme based on the original Green tensors will be used. The program is written using Matlab and all the simulations are run on a laptop computer under Windows 10 OS with Intel Core i5-10310U CPU 1.70GHz RAM 16 GB.

In the first problem, the gas is inside a parallel slit pore of width 2a = 2/3 and the characteristic length is fixed κ g /ωρ g c g = 1 (or δ = 1/ √ i). After 39 iterations which takes 1.00 s, the algorithm converges, yielding results given in Fig 3 . From the real temperature field, it is zero outside the gas phase and increases around the center due to the effect of the source term. The temperature profile cutting across the unit cell confirms this trend and shows that both the real and imaginary parts of the temperature match perfectly the analytical solutions given in the Appendix.

In the second problem, the gas is contained inside the cylindrical pore with radius a = 1/3. The characteristic length parameters are varied κ g /ωρ g c g = 1, 0.02, 0.01. The algorithm stops after 45, 78, 104 iterations (or 1.02s, 1.65s, 2.20s computation time) and results are plotted in Fig. 4. Again, the temperature is zero in the solid phase and the solution in the gaseous phase agrees with the analytical solutions. It is noted that when κ g /ωρ g c g and δ decreases, the behavior corresponds to a lower diffusivity/conductivity and the temperature profile tends to decrease and flatten out. Theoretically, the temperature field tends to a vanishing small field θ → 0 when δ → 0. When δ increases, the conductive behavior becomes increasingly dominant. Numerical evidence shows that when κ g /ωρ g c g ≥ 1, the temperature profile is stabilized and less sensitive to the parameter κ g /ωρ g c g . Figure 4: The case of circle pore with perfect solid-gas interface. The real part of temperature field ℜ(κ g θ/iωP ) with κ g /ωρ g c g = 1 (top) and the variation of the real part of temperature on axis x 2 different value of κ g /ωρ g c g (bottom).

Next we shall consider a problem on a fiber material containing a random setting of cylindrical fibers. The sample section containing 15 random nonoverlapping circles of radius 0.06 is generated by the Random Sequential Addition (RSA). The cylinder represents the solid phase and the pore corresponds to the space between the cylinders. We are interested by the convergence behavior of the problem, i.e the variation of errors ε je and ε θ with respect to the number of iterations. From Fig. 5, the error norms exhibit a common two regimes behavior: a steep descent followed by a slower descent. This two regimes behaviour is also found in the recent work on the convergence of FFT based iterative schemes [START_REF] Moulinec | Convergence of iterative methods based on Neumann series for composite materials: theory and practice[END_REF]. It is noted that our error norm exhibits some curvature rather than straight lines in [START_REF] Moulinec | Convergence of iterative methods based on Neumann series for composite materials: theory and practice[END_REF]. It is suggested that implementation of temperature constraints has some effects on the variation trend of the error norm. While the slope of the slow descent is relatively similar for all error norms, the slope and the extent of the steep descent depend on κ s . Using high values of κ s accelerates the convergence and leads to a small error value earlier and before starting the slow regime.

permeability problem with temperature jump boundary

In this subsection, the performance of the algorithm for imperfect boundary condition will be studied. Let us study again the first two examples, parallel pore and circle pore with the same parameters as in the previous section and introduce the interface coefficient l = δ = 1/ √ i (or κ int = κ g /δ) at the solid gas boundary. The global tolerance is the same as the previous case ε = 0.001. We note that for normal homogenization problems of heterogeneous materials, to achieve the best convergence rate, the reference conductivity κ 0 should be taken as κ 0 = √ κ min κ max and the two positive values κ max , κ min are the upper and lower bounds of the local conductivity.

Under the special problem with presence of the interfaces, the void phase and the external source at the interface, the optimal reference conductivity κ 0 , is generally unknown and must be studied case by case. In the slit pore problems, it is found that the FFT method has a good performance. Specifically, with κ 0 = κ g , we need 4 iterations (0.15 s) to reach the global accuracy ε and the results compare well with the analytical solution.

Both imaginary and real parts produce the jump value at the interface and the variation in the pore is also the same (see Fig. 6). In the figure, the temperature in the solid phase obtained by FFT method is not zero because we are interested only in the gas temperature. The solid conductivity is set to zero (κ s = 0) and the heat exchange between the solid and the gas is accounted for by the source term. It is noted that for slit pore, the real interface nearly coincides with the pixel based boundary at the resolution N = 128.

For the cylindrical pore problem (see Fig. 7), we observe a significant difference between the scheme without curve correction (42) and the one using pixel based approximation (43). So, the comparison ls made also with curve correction in the source term. Without the correction, the resulting source term which is proportional to the measure of the boundary can be over or underevaluated. This results in a deviation from the analytical solution. When taking account for the length of the circle crossing the cell, the agreement is excellent. The convergence is also fast as with κ 0 = κ g , the algorithm stops after 15 iterations (or 0.35 s)

It can be concluded that to study interface effects, the evaluation of the boundary flux is important and requires an accurate description of the boundary length in the cells.

Next the problem of fiber material with varying slip length l from 10 -3 δ to 10 5 δ is considered. To ensure the convergence, the reference conductivity κ 0 must also vary, especially for small values of slip length l, near an ideal interface regime. It requires an increasing value of κ 0 in these and the is slower. for l = 0.001δ, we need κ 0 = 1000κ g and 11707 iterations (or 250 s). For moderate and high slip l ≥ 0.1δ, the method based on penalty technique is very fast, with a number of iterations ranging from 10 to 200. Looking at the temperature field in Fig. 8, a significant effect of the slip length on the solution can be observed. When the slip length increases, there is no resistance at the interface and the solid phase tends to facilitate the heat flow in the gas. Consequently, the temperature field in the gas is much more homogeneous. At the limit case, the homogeneous temperature field is expected to satisfy (2) and the temperature field can be obtained

θ = - iωP κ g δ 2 (60)
and the thermal permeability Θ according to ( 9) With δ = 1/ √ i, these theoretical limits agree with the numerical results in Fig. 9

Θ(δ(ω)) = -δ 2 . ( 61
)

Comparison of the results obtained using static or frequency dependent tensor schemes

Up to now , all results have been obtained by using the "static" Green tensor. So, in this subsection, the comparison is made between the results obtained from both Green tensors. A firt result is that the frequency dependent tensor schemes P ′ , Q ′ , R ′ , H ′ yield the same solution at convergence as the original scheme for all the problems considered in the previous subsection. In terms of convergence rates, they are obviously the same when the reference heat capacity vanishes ρ 0 c 0 = 0. To achieve the comparison, the new scheme is used by varying the reference parameter ρ 0 c 0 , which can be an For the problem of slit pore (see Table 3), a behavior similar to the one of cylindrical pore is observed. There is some small difference in the fact that the scheme is stable for a larger range of ρ 0 c 0 ρgcg , especially when ρ 0 c 0 ρgcg is real or not purely imaginary. The scheme converges relatively well up to | ρ 0 c 0 ρgcg | = 100 -1000, compared to | ρ 0 c 0 ρgcg | = 0.6 -0.7 for the cylindrical pore. This is a special problem where the heat transfer is along one direction only. It suggests that the particularity of the solution renders the spectral property of the scheme less sensitive to ρ 0 c 0 . From these two examples, we find that the dynamic Green tensor, though being a priori more "adapted" to the dynamic equation, does not improve in any way the speed or convergence, compared to the static Green tensor. The best convergence being observed when the heat capacity is set to zero in the dynamic Green tensor.

ρ 0 c 0 ρgcg N iter ρ 0 c 0 ρgcg N iter ρ 0 c 0 ρgcg N iter -

Concluding remarks

In the present paper, the Fourier Transform method is developed to solve the frequency dependent heat transfer problem. The main ingredient of the method is the LS equation for heat transfer problem based on discrete Green operators. The source term appearing in the equation is used to capture the harmonic fluctuation of pressure and to enforce the constraints at the boundary. Both ideal and temperature jump boundary conditions can also be modelled by the method. From the applications to different microstructures, the method proves to be fast and to yield accurate results.

A generalized scheme based on a dynamic Green tensors is also proposed. Analogous to elastodynamics, those frequency dependent tensors are derived by using a reference complex heat capacity ρ 0 c 0 . While being more flexible, the scheme converges equally well at small value of ρ 0 c 0 . Equivalently, the optimal convergence of the generalized scheme is recovered in the particular case where the dynamic Green tensor coincides with the static Green tensor.

The use of source terms was shown in a previous work to be very useful when solving boundary value problems using FFT formulation [START_REF] To | Fourier transform approach to non periodic boundary value problems in porous conductive media[END_REF]. This work shows that the introduction of these source terms can be also of interest when the sources are frequency-dependent as in the case of dynamic heat conduction.

Appendix A. Analytical solutions for 2D heat transfer problems with temperature jump

For parallel pore of width 2a limited by the inequality -a < x 2 < a (see Fig. A10a), the solution θ(x 2 ) in the gas phase must satisfy

θ ′′ (x 2 ) - θ(x 2 ) δ 2 + iωP κ g = 0. (A.1)
The general solution θ = P ρ g c g + Ae x 2 /δ + Be -x 2 /δ (A.2) e x 2 /δ + e -x 2 /δ (e a/δ + e -a/δ ) + l/δ(e a/δ -e -a/δ ) .

(A.4)

In the case of perfect interface l/δ = 0, the classical solution is recovered θ(x 2 ) = P ρ g c g 1 -e x 2 /δ + e -x 2 /δ e a/δ + e -a/δ (A.5)

where temperature equal to 0 at x 2 = ±a.

For circular pore of radius a (see Fig. A10b), the governing solution in the cylindrical coordinate for θ(r) in the gas phase becomes θ ′′ (r) + θ ′ (r) r -θ(r) δ 2 + iωP κ g = 0 (A.6)

Figure 1 :

 1 Figure 1: Sketch of two phase material with high contrast conductivity ratio κ g /κ s ≪ 1.

Figure 2 :

 2 Figure 2: Boundary Γ crossing cells 1, 2, .., M .

Figure 3 :

 3 Figure 3: The case of parallel pore with perfect solid-gas interface. The real part of temperature field ℜ(κ g θ/iωP ) (top) and the variation of the real and imaginary part of temperature along axis x 2 (bottom).

Figure 6 :

 6 Figure6: The case of parallel pore with temperature jump at solid-gas interface. The variation of the real and imaginary part of temperature on axis x 2 for.

Figure 8 :

 8 Figure 8: The case of fiber material with temperature jump at solid-gas interface. The radius of the fibers is 0.06 and the porosity is ν = 0.865. The real part of temperature field ℜ(κ g θ/iωP ) at l = 0.1δ (top) and l = δ (bottom).

Figure 9 :

 9 Figure9: The real and imaginary part of effective thermal permeability Θ as function of thermal slip length l/δ. The complex length is fixed at δ = 1/ √ i.

Figure

  Figure A.10: Microstructures of two examples.

Table 1 :

 1 Different iteration schemes with unknown and operators associated to Lippmann-Schwinger equations with source term

Table 3 :

 3 Number of iteration at convergence as a function of reference heat capacity ρ 0 c 0 for slit pore problem. As for computation time, it varies from 1.0 s for 39 iterations to 2.50s for 118 iterations.

Error ( s =10)

Error j,e ( =100)

Error ( s =100)

Error j,e ( s =1000)

Error ( s =1000)

Figure 5: The case of fiber material with perfect solid-gas interface. The radius of the fibers is 0.06 and the porosity is ν = 0.8221. The real part of temperature field ℜ(κ g θ/iωP ) (top) and the variation error norms with respect to the number of iteration (bottom). arbitrary complex number. Let us pose ρ 0 c 0 /ρ g c g = re iφ with r being a real number and φ = 0, π/4, π/2. The problems of slit pore and cylindrical pore with perfect boundary conditions and parameters κ s = 100, κ g /ωρ g c g = 1 are considered as examples.

From Table 2 for cylindrical pore problem, it is found that the convergence is sensitive to the reference parameter. For all cases under consideration, the best convergence is achieved when ρ 0 c 0 = 0 or | ρ 0 c 0 ρgcg | is small. When ρ 0 c 0 ρgcg reaches a certain threshold value, the error norm decreases very slowly or even increases (divergence). The convergence behavior is symmetric if ρ 0 c 0 ρgcg is real or in the form ρ 0 c 0 ρgcg = re iπ/4 . The number of iterations remains relatively low and then changes drastically at | ρ 0 c 0 ρgcg | = 0.6 -0.7. The behavior is different when ρ 0 c 0 ρgcg is purely imaginary. The scheme still converges well at much higher value | ρ 0 c 0 ρgcg | and diverges when ρ 0 c 0 ρgcg = 100i. It is suggested that this behavior can be linked to the singularity of the new Green tensors. It is noted that the term iωρ 0 c 0 /κ 0 can be negative and the denominator |k| 2 + iωρ 0 c 0 /κ 0 of the new Green tensors can be singular for some values of k.

The general solution has the form

where I k is the modified Bessel function of first kind and of kth order. Under the boundary conditions at the gas-solid interface is given by

The exact solution of the above problem is θ(r) = P ρ g c g 1 -I 0 (r/δ) I 0 (a/δ) + l/δI 1 (a/δ) . (A.9)

In the case of perfect interface l/δ = 0, the classical solution is recovered θ(r) = P ρ g c g 1 -I 0 (r/δ) I 0 (a/δ) (A.10)

where temperature equal to 0 at r = a.


Denis Merigoux

Experience report: implementing a real-world, medium-sized program derived from a legislative specification

Keywords: legal expert systems, catala, hosuing benefits, legal formalization

INTRODUCTION

Implementing computer programs from legislative specifications has been a longstanding academic endeavor [START_REF] Layman E Allen | Symbolic logic: A razor-edged tool for drafting and interpreting legal documents[END_REF][START_REF] Trevor | Logic programming for large scale applications in law: A formalisation of supplementary benefit legislation[END_REF][START_REF] Dewey | Logical method and law[END_REF][START_REF] Knaebel | Une galère pas possible » : quand la Caf refuse de prendre en compte la résidence alternée[END_REF][START_REF] Mccarty | A Language for Legal Discourse I. Basic Features[END_REF][START_REF] Thorne | Some arguments about legal arguments[END_REF][START_REF] Sergot | The British Nationality Act As a Logic Program[END_REF]]. However, few of these experiments has so far sought to replicate a computer program that is already in production in public administration; most of them have created new programs or formalizations that cover sections of the law not previously automatically enforced by public administration. A somewhat dated but accurate state of the art of the use of legal expert systems in government agencies was compiled by [START_REF] Oskamp | AI in Law Practice? So far, not much[END_REF].

Building on the work around the Catala domain-specific language [START_REF] Huttner | Catala: Moving Towards the Future of Legal Expert Systems[END_REF]Merigoux et al. 2021a], we have chosen to engage in an exercise of replication of the existing IT system that computes the French housing benefits (and various other benefits) within the CNAF administrative agency, CRISTAL. More general context and main non-technical findings of this replication exercise can be found in [START_REF] Merigoux | Rules, Computation and Politics: Scrutinizing Unnoticed Programming Choices in French Housing Benefits[END_REF]; this presentation will focus on the challenges and lessons learned about the programming act itself, in an effort of consolidation of knowledge for this line of research.

Housing benefits in France are distributed to approximately 5 million households, for a yearly total amount of about 18 billion euros. Each household must declare its family and housing situation to the CNAF or CCMSA administrative agencies, that are tasked with computing and distributing the correct amount of housing benefits entitled to this household. The rules for computing the amount of benefit entitled from the family and housing situation are scattered in various legislative and regulatory documents, ranging from the Construction and Housing Code to the executive order of September 27 th , 2019 relative to the computation of personal housing benefits.

As permitted by French law, we asked the source code of CRISTAL in 2020 and obtained a CD-ROM containing a snapshot of COBOL sources for approximately 6 million lines of code, similarly to the previous experiment of Berne [2018]. These COBOL files, lacking any documentation about how to compile or run them, appear to be generated by a Computer-Aided Software Engineering (CASE) tool, whose high level sources have not been published. In other terms, what the CNAF published is merely the low-level source code, not meant to be read or edited by humans, while the real high-level source edited through a CASE tool is probably still hidden behind the CNAF curtains. Moreover, we received in November 2021 during subsequent discussion with the CNAF another excerpt from the source of the housing benefits computation. This excerpt was written using Oracle's proprietary Oracle Intelligent Advisor (OIA) programming language, as it appears that the CNAF had recently undertaken a rewrite of the housing benefits computation algorithm using this new technology, as a part of a grand scheme to modernize the CRISTAL system. We asked for the rest of the source code written in OIA but got no answer.

Author's address: Denis Merigoux, denis.merigoux@inria.fr, Inria, Paris, France.

Being unable to exploit the existing official codebase, we decided to replicate it. The motivation for tackling this particular piece of legislation comes from frequent media occurrences of dysfunctions in the CRISTAL system [START_REF] Knaebel | Une galère pas possible » : quand la Caf refuse de prendre en compte la résidence alternée[END_REF][START_REF] Monnet | Qui ne demande rien n'a rien" : de mauvaise volonté, la CAF vous prive peut-être injustement de ces prestations[END_REF][START_REF] Zerouala | La réforme des APL vire au cauchemar pour les allocataires et ses agents[END_REF]. We believe that the replication of this system using the Catala language and methodology can help improve the transparency and correctness of the existing system, by comparing its output with the existing CRISTAL system. We have already collected our recommendation to the CNAF regulators in an administrative report [START_REF] Merigoux | Observations sur le calcul des aides au logement[END_REF]].

The contributions of this presentation are the following: a detailed account of the programming methodology used, a summary of the legal difficulties encountered, and a deep focus on the semantic and language constructs needed to correctly express the computational content of the French housing benefit legislation.

PROGRAMMING METHODOLOGY

In absence of an exploitable codebase, we decided to start the programming from scratch and from the official legal sources describing the computation of the housing benefits, following the Catala methodology mentioned in [START_REF] Huttner | Catala: Moving Towards the Future of Legal Expert Systems[END_REF]. The main authors of the program include a programmer very proficient in the Catala language (A) and a lawyer initially specialized in intellectual property law (B). Another programer (C) switched with A for a small part of the development. None of the three had any prior experience of the French housing benefits computation.

The development effort spanned 9 months from December 2021 to August 2022, but the developers were only working part-time on the implementation during this period, as administrative and publication-related tasks take up a significant amount of time in a research lab. The time spent in programming and programming-related tasks has been self-reported in order to produce accurate estimates of the effort involved in the project. These estimates do not include the research and development time invested in the building of the Catala compiler and tooling; the goal is to measure the marginal cost of implementing a real-world program derived from legislative specifications.

The first step of the development was to gather all the legislative and regulatory sources describing the computation of the French housing benefits. This was done by B in approximately 20 hours of work, resulting in a 106-pages document compiling the (almost, see §3) complete specification.

Then, the bulk of the programming was performed in 22 sessions of pair programming totalling 50 hours. In each session, A or C was paired with B. During these sessions, intense cross-disciplinary discussion was necessary to build a shared knowledge of the effects of the law and the Catala program. But we can say in retrospect that the majority of the time was spent on actually deciphering and understanding the text of the law, as the drafting style is very compact and sometimes arcane. The Catala code produced amounts to about 7,000 lines of code. The ratio of lines of code per hour of programming is very high because of several factors. First, the code contains a lot of copy-pasting that mirrors the duplications in the law. Second, a lot of boilerplate code comes from the translation of huge tables of values, contained in an executive order assigning threshold amounts of benefits. The logic of the table had to be expressed in terms of conditionals and pattern matching in a very repetitive fashion. Third, programmer A being the very proficient the language, this scenario is one of maximum programming productivity; the project time is likely to go up with another programmer less accustomed with the language.

Next, approximately 50 hours were spent testing and debugging the resulting program. Some tests were crafted manually by B who invented fictional households and manually computed the housing benefits amounts by applying the rules in the legislative sources. Other test cases included expected outputs were collected from the official explainer document for the housing benefits computation from the government. These official test cases yielded some differences with our computations, so we investigated with the authors of the official government document, the DGALN/DHUP/FE4 administrative office in the Ministry responsible for Housing. These investigations around the test cases revealed a number of bugs in our implementation, as well as some factual errors in the official test cases that were consequently fixed. The total number of test cases involved in the qualification of the program is quite low, around 20. Indeed, crafting a new test case is very time consuming due to the complexity of the input space and the tedious nature of the housing benefits computations that involve many steps when done by hand. We do not believe the current level of testing to be sufficient and a doubling or tripling of the testing time would be absolutely necessary to reach a minimum level of assurance.

Additional software validation was performed using the initial features of the Catala proof platform [START_REF] Delaët | Turning Catala into a Proof Platform for the Law[END_REF], including the formal verification of the well-behaved execution of the program (there is always one rule that apply, no two exceptions overlap). Even though the proof platform is not fully implemented, with only a Z3 backend and a crude encoding that produces a lot of false positives, it found a bug. Finally, we developed a tool that scans our codebase and compares the legislative text with the official version stored on the Légifrance database. This tool helped us locate and fix a few copy-pasting errors (especially in tables) that were mirrored in the Catala code.

LEGAL CHALLENGES

As reported by [START_REF] Merigoux | Rules, Computation and Politics: Scrutinizing Unnoticed Programming Choices in French Housing Benefits[END_REF], we found a number of factual errors in the text of the executive orders describing the computation parameters for the housing benefits. The errors included two typos in numbers contained in a table of parameters directly determining the maximum amount of benefits entitled in a corner case ("59.47" instead of "459.47" and "47.99" instead of "347.99"), and more problematically a whole table of income brackets was missing for a mode of computation (although it was published in the official, non-legal explainer document from the government). Those factual errors have been fixed by the government with a new executive order amending the previous one.

Another legal challenge we encountered was the lack of clarity in the drafting style of executive orders. Indeed, as reported by the the DGALN/DHUP/FE4 administrative office, the executive orders describing the French housing benefits computation had undergone in 2019 a massive restructuring and reform, that led to a coexisting of old and new texts in the current version. The articulation between the old and new texts had not yet been perfectly consolidated at the time of our development, and we sent a list of a dozen questions to the DGALN/DHUP/FE4 office to clarify the interpretations. DGALN/DHUP/FE4 answered promptly and later requested our input for a coming "grooming" rewrite of these texts, which has not yet been enacted at submission time.

The computer replication of the French housing benefits computation, conducted with a formal and detailed-oriented methodology, led to two interesting findings at the intersection between law and computation. First, as detailed by [START_REF] Merigoux | Observations sur le calcul des aides au logement[END_REF], a computational trick in the final steps of the computation only works if the rate of the CRDS tax is below 1%. If the rate is ever changed to a value superior to 1% (currently it is 0.5%), the formula written in the executive order will result in rounding errors potentially detrimental to the beneficiaries. This finding stresses the importance of making computational invariants explicit, even in legislative texts. Second, the computation distinguishes between two exceptional cases for the housing benefits amount: the case where the housing is a room within a larger housing unit, and the case where the housing is shared with flatmates. However, ambiguity prevails when the two situations happen at the same time. Confronted with that ambiguity, DGALN/DHUP/FE4 answered: "As this situation is far from reality, it has not been dealt with". This reminds us of a delayed form of decision that clashes with the nature of computer programs [START_REF] Diver | Interpreting the Rule(s) of Code: Performance, Performativity, and Production[END_REF][START_REF] Hildebrandt | Code-driven Law: Freezing the Future and Scaling the Past[END_REF]. To be respectful of that delayed decision, the best option could be for the program to crash and instruct a human review of the case.

The last but not least of the legal challenges was the computation of the housing benefits in the case of split custody. [START_REF] Merigoux | Observations sur le calcul des aides au logement[END_REF] provides an in-depth explanation of the thorny situation on this issue: pressed by past court decisions to split the benefits between the two households of the split custody child, the administrative agencies delay the application of this measure because of a technical blocking point. This technical blocking point is in fact related to the capacity of the IT system to compute the benefits twice for the same household, with or without the children in split custody. The ability to call the computation twice requires proper computer programming abstractions, like functions or classes, that have not yet been identified as a key requirement for legal expert systems (see §4). Our Catala replication has proven flexible enough to express the exact computation of the benefits split, therefore providing a technical path to resolving the stillstand that has been ongoing since 2017. Since then, only households that file an official complaint to the administration have their benefits split by hand, the others cannot have their split custody child accounted for properly.

LANGUAGE DESIGN CHALLENGES

In this section, we will discuss the language design decisions of Catala with respect to the empirical needs of the French housing benefits computation.

Defeasability. The need for a form of defeasible logic when formalizing legislation has been identified for a long time, as summarized by [START_REF] Thorne | Some arguments about legal arguments[END_REF]]. Catala embarks a very restricted form of defeasible logic in its default calculus [Merigoux et al. 2021a], corresponding to a subset of prioritized default logic [START_REF] Brewka | Prioritizing default logic[END_REF]. We found that this restricted form of defeasible logic was sufficient to correctly express all the exceptional case structures found in the French housing benefits computation. The most complex exceptional case structure involved multiple groups of rules, each group piecewise defining the value of a variable, being arranged in a tree-like structure where a group of rules can have multiple groups of rules as exceptions to it, and the exceptions can also recursively have exceptions. We found the lean semantics and runtime for dealing with exceptions in Catala to be sufficient for our needs.

Precision. The computation of the French housing benefits involves amounts of money being multiplied and divided multiple times. Particularly, the computation of the "rent equivalence" involves up to 6 income brackets and a division by twelve that is sufficiently complex to raise the usual problems of computing precision stemming from the binary representation of numbers in the machine. Learning from the mistakes of the French tax authority that chose floating-point numbers to represent all variables, including amounts of money [Merigoux et al. 2021b], Catala features a rich typology of values where amounts of money, integers and decimal numbers each have different representations. In particular, the current Catala semantics represents money as an integer number of cents; multiplying a money by a decimal yields a new amount of money rounded to the nearest cent. Hence, in an expression like (𝑚 × 𝑎) × 𝑏 where 𝑚 is an amount of money and 𝑎 and 𝑏 two decimals, two roundings to the nearest cent happen when computing with Catala. In the detailing of the official tests cases for French housing benefits as published by DGALN/DHUP/FE4, all intermediate amounts of money are also rounded to the nearest cent. But when confronted with the question of computing precision, DGALN/DHUP/FE4 later clarified that the rounding to the nearest cent was merely an "explanation artefact", and that large formulae in the computation should only round at the end and not in intermediate values. We had to change our Catala code (switching from the money type to decimal) to account for that behavior, that could change by up to 1 euro the amount of benefits effectively received. More generally, the question of computing precision for taxes and social benefits computation seems to be largely under-considered by both lawyers and programmers, compared to the potentially massive cumulated effects it can have.

Modules and functions. The computation of the French housing benefits, with 7,000 lines of code, is big enough to trigger some software engineering challenges. First, we felt the need for modules representing separate compilation units to better structure our code. This limitation is being addressed with a planned module system for Catala. Second, there are six different modes of computation with four main formulae that can interact or reference each other. Because of these interactions, being able to encapsulate these modes of computation and formulae into functions of our language was paramount. To that end, Catala features the concept of "scope", that corresponds to a function with some defeasible capabilities builtin. The current Catala syntax only allows a static declaration of the scope call graph, which proved to be a limitation since we want certain sub-scopes to be called only if a condition is filled. The addition of a reified scope call syntax that lifts this limitation is underway. The functional nature of Catala with immutable values lead to a certain amount of boilerplate "plumbing" code to fill the inputs of sub-scopes with values that are themselves inputs of the parent scope. However, this explicitation of the input values of the scope has proved comfortable when dealing with subtle counter-factuals in the law that mandate a certain computation to be made again with slightly different arguments. The usefulness of functions for formalizing legislation is somewhat under-represented in the literature, and our findings here are consistent with previous studies on the French tax code implementation [Merigoux et al. 2021b].

Higher-order. Also somewhat under-represented in the literature, the usefulness of higher-order functions has been empirically validated in the French housing benefits computation. Used in many places, higher-order functions felt natural to the lawyer that reviewed the code. In Catala, functions only have one argument so we don't run into partial application or curryfication problems. Even though a higher-order function usually can be replaced by a refactoring, they were critical in a high-stakes situation. For splitting benefits in the case of split custody, what has to be split is a "gross" amount of benefits, on which several post-processing steps have to be applied to reach the final amount. These steps depend on the mode of computation used (6 in total), but we didn't want to have 6 copies of the splitting operation that is inserted before the post-processing steps. By reifying the post-processing steps with a higher-order function, we were able to avoid code duplication and correctly express the computation.

CONCLUSION

This case study on the French housing benefits has been conducted not only with the objective of producing a piece of software, but also to interrogate the ways in which this particular piece of law was written and enforced by the government. The numerous findings it yielded can inform both lawyers and computer scientists working at the intersection of their fields.

This replication provides a starting point for a potential replacement as the production system for computing the French housing benefits. Even though the choice of replacing a production system is consequential and involves a lot of parameters, we believe that our replication and the underlying Catala technology have the potential of being production-ready for government systems within a few years.

ACKNOWLEDGMENTS

Many thanks to Lilya Slimani who has been my pair programming partner for the development of the French housing benefits codebase, and to Alain Delaët for the occasional help. Thanks to the DGALN/DHUP/FE4 administrative office for their very reactive feedback and cooperation, which helped us debug our program and somewhat improve the redaction of the executive orders describing the computation of the French housing benefits.