# 3D modeling of an arc in interaction with a composite panel

François Pechereau, Fabien Tholin, Julien Vanharen, Guillaume Puigt, Frédéric Alauzet, Philippe Lalande

► To cite this version:

François Pechereau, Fabien Tholin, Julien Vanharen, Guillaume Puigt, Frédéric Alauzet, et al. 3D modeling of an arc in interaction with a composite panel. Aerospace Europ Conference - 10th EU-CASS - 9th CEAS Conference, Jul 2023, Lausanne, Switzerland. 10.13009/EUCASS2023-701. hal-04401291

HAL Id: hal-04401291

https://inria.hal.science/hal-04401291

Submitted on 17 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# 3D modeling of an arc in interaction with a composite panel

François Pechereau, Fabien Tholin, Julien Vanharen, Guillaume Puigt, Frédéric Alauzet, Philippe Lalande

► To cite this version:

François Pechereau, Fabien Tholin, Julien Vanharen, Guillaume Puigt, Frédéric Alauzet, et al. 3D modeling of an arc in interaction with a composite panel. Aerospace Europ Conference - 10th EU-CASS - 9th CEAS Conference, Jul 2023, Lausanne, Switzerland. 10.13009/EUCASS2023-701. hal-04401291

HAL Id: hal-04401291

https://inria.hal.science/hal-04401291v1

Submitted on 17 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 3D modeling of an arc in interaction with a composite panel

François Pechereau*, Fabien Tholin*, Julien Vanharen†, Guillaume Puigt°, Frédéric Alauzet† and Philippe Lalande*

*DPHY, ONERA, Université Paris Saclay

F-91123 Palaiseau-France

†GammaO, Equipe Inria/Onera, Centre Inria de Saclay, Université Paris-Saclay

F-91120 Palaiseau-France

°GammaO, Equipe Inria/Onera, ONERA/DMPE, Université de Toulouse

F-31055 Toulouse-France

francois.pechereau@onera.fr · fabien.tholin@onera.fr

*Corresponding author

## Abstract

During the lightning strike on an aircraft and particularly in the pulsed arc phase, the current can reach values in the range of 200 kA during microsecond timescale. The direct effect of such arcs on structures can be quite severe depending on the material being used. This effect is even more troublesome on new composite materials used in modern airplane such as such as multilayer CFRP (Carbon Fiber Reinforced Polymer). Indeed, due to the highly anisotropic properties of CFRP composites, the electrical conductivity of composites is not scalar as in metallic materials but is generally tensorial. In a multilayer CFRP, each layer exhibits an electrical conductivity at least two order of magnitude higher in the direction of the fibers than in the direction perpendicular to it. Moreover, the different layers have different orientations of the fibers, and the transverse conductivity between them may be three or four orders of magnitude lower than the conductivity of the fibers. For the numerical study of an arc in the pulsed phase interacting with CFRP composite materials then requires to perform 3D simulations. After the presentation of the MHD code Taranis and the anisotropic mesh refinement code geflo.a, the new Taranis/geflo.a framework will be used to carry out a 3D numerical study showing the effect of the tensorial electrical conductivity of a composite on the dynamics of an arc in the pulsed phase. The simulation results will then be compared with experimental results.

## 1. Symbols

$\vec{j}$: Current density vector [A.m⁻²],

$\sigma$: Electrical conductivity [S.m⁻¹],

$\phi$: Electrical potential [V],

$\vec{E}$: Electric field [V.m⁻¹],

$\vec{B}$: Magnetic field [T],

$\vec{A}$: Magnetic vector potential [T.m],

$\rho$: Air density [kg.m⁻³],

$\vec{u}$: velocity vector [m.s⁻¹],

P: Pressure [Pa],

E: Total specific energy [J.kg⁻¹],

$S_{rad}$: Radiation energy source term [W.m⁻³],

$\kappa_v$: Absorption coefficient [m⁻¹] of air at frequency $v$ [Hz],

$G_v$: Spectral incident flux [W.m⁻².Hz⁻¹] at frequency $v$,

$L_v^0$ (T): Black body spectral emittance [W.m⁻²] at frequency $v$ and temperature T [K].

## 2. Introduction

During the pulsed arc phase of a lightning strike to an aircraft, the current may reach values as high as 200 kA during microsecond timescales (Fisher and Plumer, 1977). In previous studies (Tholin et al., 2015) and (Tholin et al., 2017), a
# 3D modeling of an arc in interaction with a composite panel

François Pechereau*, Fabien Tholin*, Julien Vanharen†, Guillaume Puigt°, Frédéric Alauzet† and Philippe Lalande*

*DPHY, ONERA, Université Paris Saclay

F-91123 Palaiseau-France

†GammaO, Equipe Inria/Onera, Centre Inria de Saclay, Université Paris-Saclay

F-91120 Palaiseau-France

°GammaO, Equipe Inria/Onera, ONERA/DMPE, Université de Toulouse

F-31055 Toulouse-France

francois.pechereau@onera.fr · fabien.tholin@onera.fr

*Corresponding author

## Abstract

During the lightning strike on an aircraft and particularly in the pulsed arc phase, the current can reach values in the range of 200 kA during microsecond timescale. The direct effect of such arcs on structures can be quite severe depending on the material being used. This effect is even more troublesome on new composite materials used in modern airplane such as such as multilayer CFRP (Carbon Fiber Reinforced Polymer). Indeed, due to the highly anisotropic properties of CFRP composites, the electrical conductivity of composites is not scalar as in metallic materials but is generally tensorial. In a multilayer CFRP, each layer exhibits an electrical conductivity at least two order of magnitude higher in the direction of the fibers than in the direction perpendicular to it. Moreover, the different layers have different orientations of the fibers, and the transverse conductivity between them may be three or four orders of magnitude lower than the conductivity of the fibers. For the numerical study of an arc in the pulsed phase interacting with CFRP composite materials then requires to perform 3D simulations. After the presentation of the MHD code Taranis and the anisotropic mesh refinement code geflo.a, the new Taranis/geflo.a framework will be used to carry out a 3D numerical study showing the effect of the tensorial electrical conductivity of a composite on the dynamics of an arc in the pulsed phase. The simulation results will then be compared with experimental results.

## 1. Symbols

$\vec{j}$: Current density vector [A.m$^{-2}$],

$\sigma$: Electrical conductivity [S.m$^{-1}$],

$\phi$: Electrical potential [V],

$\vec{E}$: Electric field [V.m$^{-1}$],

$\vec{B}$: Magnetic field [T],

$\vec{A}$: Magnetic vector potential [T.m],

$\rho$: Air density [kg.m$^{-3}$],

$\vec{u}$: velocity vector [m.s$^{-1}$],

P: Pressure [Pa],

E: Total specific energy [J.kg$^{-1}$],

$S_{rad}$: Radiation energy source term [W.m$^{-3}$],

$\kappa$: Absorption coefficient [m$^{-1}$] of air at frequency $v$ [Hz],

$G_v$: Spectral incident flux [W.m$^{-2}$.Hz$^{-1}$] at frequency $v$,

$L_v^0$ (T): Black body spectral emittance [W.m$^{-2}$] at frequency $v$ and temperature T [K].

## 2. Introduction

During the pulsed arc phase of a lightning strike to an aircraft, the current may reach values as high as 200 kA during microsecond timescales (Fisher and Plumer, 1977). In previous studies (Tholin et al., 2015) and (Tholin et al., 2017), a
Magneto-Hydro-Dynamics (MHD) model was used to get insight into the fundamental mechanisms of the interaction between a lightning arc and aerospace materials. However, the validity range of the databases for plasma properties and the numerical schemes implemented could not sustain such high current values and rise-times that may result in temperatures higher than 35 000 K and pressure higher than 100 bar. For this purpose, a new MHD code called Taranis, presented in (Pechereau et al., 2019), has been developed that do not suffer those limitations for the simulation of the pulsed arc phase of lightning stroke. Taranis solves the electromagnetic fields under static or quasi-static approximations, the radiative transfer equation and the compressible Euler or Navier-Stokes equations under Local Thermodynamic Equilibrium (LTE). Under the same framework, it is able to compute accurately the current density distribution in the arc and in the conductive materials as well as the feedback of the Joule effect, the Laplace forces and the radiative emission and absorption in the air plasma. In (Pechereau et al., 2019), validation of the MHD model in 2D axisymmetric geometry was presented. But in order to fully compare the simulation of an arc in a pulsed phase interacting with composite material, a 3D approach is mandatory due to the tensorial nature of the electrical conductivity of the composite material. In order to tackle this task, the MHD code Taranis developed at ONERA has been coupled with the anisotropic mesh refinement code geflo.a developed at INRIA's team Gamma. In this article, the fruit of this collaboration between ONERA and INRIA is presented where for the first time an arc in pulsed phase was simulated using anisotropic mesh refinement.

After the presentation of the code Taranis and the MHD model, a description of the geflo.a code is given. Starting from experimental results from the thesis (Martins, 2016), simulation results showing the dynamics of an arc in a pulsed phase in interaction with a composite material is shown and explained.

## 3. Presentation of the workflow Taranis/geflo.a

### 3.1 Taranis code

Taranis is a finite-volume compressible resistive MHD code that is able to compute the complex motion of a compressible and conductive fluid such as air plasmas under the influence of electromagnetic fields and radiative fluxes. This physical process is described by the coupled system of equations 1.

$$ \frac{\partial \rho}{\partial t} + \vec{\nabla} \cdot (\rho \vec{u}) = 0, \\ \frac{\partial (\rho \vec{u})}{\partial t} + \vec{\nabla} \cdot [\rho \vec{u} \otimes \vec{u}] = -\vec{\nabla} \vec{p} + \vec{j} \times \vec{B}, \qquad (1) $$

$$ \frac{\partial(\rho E)}{\partial t} + \vec{\nabla} \cdot (\rho \vec{u} H) = \vec{j} \cdot \vec{E} + S_{rad} \\ \vec{\nabla} \cdot \vec{j} = \vec{\nabla} \cdot \sigma \vec{E} = -\vec{\nabla} \cdot \sigma \vec{\nabla} \phi = 0, \\ \vec{\nabla} \times \vec{B} = \vec{\nabla} \times \vec{\nabla} \times \vec{A} = \mu_0 \vec{j} \qquad (2) $$

$$ \vec{\nabla} \cdot (\kappa_v^{-1} \vec{\nabla} G_v) = 3\kappa_v(G_v - 4\pi L_v^0(T)) \qquad (3) $$

In a first step, Taranis solves here the compressible Euler to compute the fluid dynamics, the pressure P and the temperature T fields. Two electromagnetic source terms are accounted for: The Laplace forces $\vec{j} \times \vec{B}$ (N.m$^{-3}$) in the momentum conservation equation, and the Joule effect $\vec{j} \cdot \vec{E}$ (W.m$^{-3}$) in the energy conservation equation. In a second step, with the updated temperature and pressure fields, Taranis computes the electrical conductivity $\sigma$ needed by the current conservation equation 2. Then this elliptic equation is solved thanks to a linear solver and the electrical potential $\phi$, the current density $\vec{j}$ and the electric field $\vec{E}$ are computed. The magnetic field $\vec{B}$ is then computed by solving the Maxwell-Ampere equation reformulated in magnetic vector potentials resulting in three linear systems to solve. Finally, the new electromagnetic source-terms are updated for the next iteration. In the Taranis code, a conjugate gradient method preconditioned with a SOR method from the PETSc library is used for the linear solver (Balay et al., 1997), (Balay et al., 2023a) and (Balay et al., 2023b).

Due to the extreme temperatures ($T > 35000K$) and pressures ($P > 100$ bar) occurring in the early stages of a pulsed arc phase, the Radiative Transfer Equation (LTE) is solved in Taranis to take into account the redistribution of the energy due to the emission and the absorption of radiations in the air plasma. A 14 spectral bands P1 approximation is
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

Magneto-Hydro-Dynamics (MHD) model was used to get insight into the fundamental mechanisms of the interaction between a lightning arc and aerospace materials. However, the validity range of the databases for plasma properties and the numerical schemes implemented could not sustain such high current values and rise-times that may result in temperatures higher than 35 000 K and pressure higher than 100 bar. For this purpose, a new MHD code called Taranis, presented in (Pechereau et al., 2019), has been developed that do not suffer those limitations for the simulation of the pulsed arc phase of lightning stroke. Taranis solves the electromagnetic fields under static or quasi-static approximations, the radiative transfer equation and the compressible Euler or Navier-Stokes equations under Local Thermodynamic Equilibrium (LTE). Under the same framework, it is able to compute accurately the current density distribution in the arc and in the conductive materials as well as the feedback of the Joule effect, the Laplace forces and the radiative emission and absorption in the air plasma. In (Pechereau et al., 2019), validation of the MHD model in 2D axisymmetric geometry was presented. But in order to fully compare the simulation of an arc in a pulsed phase interacting with composite material, a 3D approach is mandatory due to the tensorial nature of the electrical conductivity of the composite material. In order to tackle this task, the MHD code Taranis developed at ONERA has been coupled with the anisotropic mesh refinement code geflo.a developed at INRIA's team Gamma. In this article, the fruit of this collaboration between ONERA and INRIA is presented where for the first time an arc in pulsed phase was simulated using anisotropic mesh refinement.

After the presentation of the code Taranis and the MHD model, a description of the geflo.a code is given. Starting from experimental results from the thesis (Martins, 2016), simulation results showing the dynamics of an arc in a pulsed phase in interaction with a composite material is shown and explained.

3. Presentation of the workflow Taranis/eflo.a

3.1 Taranis code

Taranis is a finite-volume compressible resistive MHD code that is able to compute the complex motion of a compressible and conductive fluid such as air plasmas under the influence of electromagnetic fields and radiative fluxes. This physical process is described by the coupled system of equations 1.

$$
\begin{align}
\frac{\partial \rho}{\partial t} + \vec{\nabla} \cdot (\rho \vec{u}) &= 0, \nonumber \\
\frac{\partial (\rho \vec{u})}{\partial t} + \vec{\nabla} \cdot [\rho \vec{u} \otimes \vec{u}] &= -\vec{\nabla} \vec{p} + \vec{j} \times \vec{B}, \tag{1} \\
\frac{\partial (\rho E)}{\partial t} + \vec{\nabla} \cdot (\rho \vec{u} H) &= \vec{j} \cdot \vec{E} + S_{\text{rad}} \nonumber
\end{align}
$$

$$
\begin{align}
\vec{\nabla} \cdot \vec{j} &= \vec{\nabla} \cdot \sigma \vec{E} = -\vec{\nabla} \cdot \sigma \vec{\nabla} \phi = 0, \tag{2} \\
\vec{\nabla} \times \vec{B} &= \vec{\nabla} \times \vec{\nabla} \times \vec{A} = \mu_0 \vec{j} \nonumber
\end{align}
$$

$$
\vec{\nabla} \cdot (\kappa_v^{-1} \vec{\nabla} G_v) = 3\kappa_v(G_v - 4\pi L_v^0(T)) \quad (3)
$$

In a first step, Taranis solves here the compressible Euler to compute the fluid dynamics, the pressure P and the
temperature T fields. Two electromagnetic source terms are accounted for: The Laplace forces $\vec{j} \times \vec{B}$ (N.m$^{-3}$) in
the momentum conservation equation, and the Joule effect $\vec{j} \cdot \vec{E}$ (W.m$^{-3}$) in the energy conservation equation. In a
second step, with the updated temperature and pressure fields, Taranis computes the electrical conductivity $\sigma$ needed
by the current conservation equation 2. Then this elliptic equation is solved thanks to a linear solver and the electrical
potential $\phi$, the current density $\vec{j}$ and the electric field $\vec{E}$ are computed. The magnetic field $\vec{B}$ is then computed by
solving the Maxwell-Ampere equation reformulated in magnetic vector potentials resulting in three linear systems to
solve. Finally, the new electromagnetic source-terms are updated for the next iteration. In the Taranis code, a conjugate
gradient method preconditionned with a SOR method from the PETSc library is used for the linear solver (Balay et al.,
1997), (Balay et al., 2023a) and (Balay et al., 2023b).

Due to the extreme temperatures (T > 35000K) and pressures (P > 100 bar) occurring in the early stages of a pulsed
arc phase, the Radiative Transfer Equation (RTE) is solved in Taranis to take into account the redistribution of the
energy due to the emission and the absorption of radiations in the air plasma. A 14 spectral bands P1 approximation is
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

used resulting in 14 elliptic equations in the form of equation 3 to solve. The 14 spectral bands are detailed in Table 1, ranging from the mid-infrared to extreme ultraviolet. For each spectral band, a Rosseland mean of the absorption coefficient $\kappa_v$ has been considered because air plasmas in lighting arcs are optically thick. In Taranis, all the plasma properties are tabulated at LTE as a function of the temperature T and the pressure P. To be able to compute a complete pulsed arc phase with applied current peaks up to 200 kA, the data-tables of radiative, transport and thermodynamic properties have been extended to the range $300 \text{ K} < T < 70 \text{ kK}$ and $0.1 < P < 1 \text{ kbar}$. The absorption coefficients were computed by the authors of (Chauveau et al., 2003), while thermodynamic and transport properties were taken from (D'angola et al., 2008).

<table><caption>Table 1: Spectral bands used in the radiative transfer model</caption><thead><tr><th>Band</th><th>Frequency [10<sup>14</sup> Hz]</th><th>Wave length [µm]</th></tr></thead><tbody><tr><td>1</td><td>0.29-0.57</td><td>10.00-5.26</td></tr><tr><td>2</td><td>0.57-2.95</td><td>5.26-1.01</td></tr><tr><td>3</td><td>2.95-3.83</td><td>1.01-0.78</td></tr><tr><td>4</td><td>3.83-7.48</td><td>0.78-0.40</td></tr><tr><td>5</td><td>7.48-9.88</td><td>0.40-0.30</td></tr><tr><td>6</td><td>9.88-14.3</td><td>0.30-0.21</td></tr><tr><td>7</td><td>14.3-20.9</td><td>0.21-0.14</td></tr><tr><td>8</td><td>20.9-22.9</td><td>0.14-0.13</td></tr><tr><td>9</td><td>22.9-29.1</td><td>0.13-0.10</td></tr><tr><td>10</td><td>29.1-44.6</td><td>0.10-0.07</td></tr><tr><td>11</td><td>44.6-59.9</td><td>0.07-0.05</td></tr><tr><td>12</td><td>59.9-69.9</td><td>0.05-0.04</td></tr><tr><td>13</td><td>69.9-79.9</td><td>0.04-0.04</td></tr><tr><td>14</td><td>79.9-89.93</td><td>0.04-0.03</td></tr></tbody></table>

In this work the current waveform used is the one corresponding to the GRIFON lightning generator used at Onera that was used in the experiments. The characteristics of the Grifon wave noted "G" as well as the normative "D" and "A" waves are found on Table 2 with the maximum current values imax (kA), the maximum current times tmax (µs), the action integrals AI (kJ ω⁻¹) and the maximum current rise-times (di/dt)max (kA µs⁻¹).

<table><caption>Table 2: Characteristics of different current waves.</caption><thead><tr><th>Wave</th><th>i<sub>max</sub> [10<sup>3</sup>A]</th><th>t<sub>max</sub> [µs]</th><th>AI [kJ ω<sup>-1</sup>]</th><th>(di/dt)<sub>max</sub> [kA µs<sup>-1</sup>]</th></tr></thead><tbody><tr><td>G</td><td>100</td><td>12.22</td><td>130</td><td>16.8</td></tr><tr><td>D</td><td>100</td><td>3.17</td><td>250</td><td>139.1</td></tr><tr><td>A</td><td>200</td><td>6.36</td><td>2000</td><td>139.1</td></tr></tbody></table>

To initiate the simulations, a hot 2 mm wide channel with a 10 000 K Gaussian temperature profile is used to ensure a non-zero conductivity at t=0 s, that is representative of the leader channel preceding the lightning stroke. It has been checked that the total energy in the initialization channel is very small compared to the total energy input during the lightning and that it has no significant influence on the results presented in this study.

For the modelling of composite materials and specifically the electrical conductivity, implementing in Taranis a tensorial electrical conductivity model was required. A modification in the equation 2 is done where the electrical conductivity is no longer a scalar variable $\sigma$ but a tensorial variable $\overline{\sigma}$. Its expression is given by equation 4.

$$ \overline{\overline{\sigma}} = \sigma_{\parallel} \cdot \overline{\overline{C}} = \sigma_{\parallel} \cdot \begin{bmatrix} \cos\theta^2 + \frac{\sigma_p}{\sigma_{\parallel}} \sin\theta^2 & \cos\theta \sin\theta \left[1 - \frac{\sigma_p}{\sigma_{\parallel}}\right] & 0 \\ \cos\theta \sin\theta \left[1 - \frac{\sigma_p}{\sigma_{\parallel}}\right] & \sin\theta^2 + \frac{\sigma_p}{\sigma_{\parallel}} \cos\theta^2 & 0 \\ 0 & 0 & \frac{\sigma_p}{\sigma_{\parallel}} \end{bmatrix} \quad (4) $$

The $\theta$ angle is formed between the x axis and the direction chosen for the composite fibers. In this work $\theta = 0$, and the x direction is chosen to be the direction of the carbon fibers and the y direction is chosen to be the direction perpendicular to the carbon fibers. The tensorial electrical conductivity depends on the value of the electrical conductivity along carbon fibers $\sigma_{\parallel}$, on the one perpendicular to the fibers $\sigma_{\perp}$ and the electrical conductivity in the direction transverse to the fibers through the material taken here in the z direction $\sigma_{\tan}$.
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

used resulting in 14 elliptic equations in the form of equation 3 to solve. The 14 spectral bands are detailed in Table 1, ranging from the mid-infrared to extreme ultraviolet. For each spectral band, a Rosseland mean of the absorption coefficient $\kappa_v$ has been considered because air plasmas in lighting arcs are optically thick. In Taranis, all the plasma properties are tabulated at LTE as a function of the temperature T and the pressure P. To be able to compute a complete pulsed arc phase with applied current peaks up to 200 kA, the data-tables of radiative, transport and thermodynamic properties have been extended to the range $300 \text{ K} < T < 70 \text{ kK}$ and $0.1 < P < 1 \text{ kbar}$. The absorption coefficients were computed by the authors of (Chauveau et al., 2003), while thermodynamic and transport properties were taken from (D'angola et al., 2008).

<table><caption>Table 1: Spectral bands used in the radiative transfer model</caption><thead><tr><th>Band</th><th>Frequency [10<sup>14</sup> Hz]</th><th>Wave length [μm]</th></tr></thead><tbody><tr><td>1</td><td>0.29-0.57</td><td>10.00-5.26</td></tr><tr><td>2</td><td>0.57-2.95</td><td>5.26-1.01</td></tr><tr><td>3</td><td>2.95-3.83</td><td>1.01-0.78</td></tr><tr><td>4</td><td>3.83-7.48</td><td>0.78-0.40</td></tr><tr><td>5</td><td>7.48-9.88</td><td>0.40-0.30</td></tr><tr><td>6</td><td>9.88-14.3</td><td>0.30-0.21</td></tr><tr><td>7</td><td>14.3-20.9</td><td>0.21-0.14</td></tr><tr><td>8</td><td>20.9-22.9</td><td>0.14-0.13</td></tr><tr><td>9</td><td>22.9-29.1</td><td>0.13-0.10</td></tr><tr><td>10</td><td>29.1-44.6</td><td>0.10-0.07</td></tr><tr><td>11</td><td>44.6-59.9</td><td>0.07-0.05</td></tr><tr><td>12</td><td>59.9-69.9</td><td>0.05-0.04</td></tr><tr><td>13</td><td>69.9-79.9</td><td>0.04-0.04</td></tr><tr><td>14</td><td>79.9-89.93</td><td>0.04-0.03</td></tr></tbody></table>

In this work the current waveform used is the one corresponding to the GRIFON lightning generator used at Onera that was used in the experiments. The characteristics of the Grifon wave noted "G" as well as the normative "D" and "A" waves are found on Table 2 with the maximum current values imax (kA), the maximum current times tmax (μs), the action integrals AI (kJ ω⁻¹) and the maximum current rise-times (di/dt)max (kA μs⁻¹).

<table><caption>Table 2: Characteristics of different current waves.</caption><thead><tr><th>Wave</th><th>i<sub>max</sub> [10<sup>3</sup>A]</th><th>t<sub>max</sub> [μs]</th><th>AI [kJ ω<sup>-1</sup>]</th><th>(di/dt)<sub>max</sub> [kJ μs<sup>-1</sup>]</th></tr></thead><tbody><tr><td>G</td><td>100</td><td>12.22</td><td>130</td><td>16.8</td></tr><tr><td>D</td><td>100</td><td>3.17</td><td>250</td><td>139.1</td></tr><tr><td>A</td><td>200</td><td>6.36</td><td>2000</td><td>139.1</td></tr></tbody></table>

To initiate the simulations, a hot 2 mm wide channel with a 10 000 K Gaussian temperature profile is used to ensure a non-zero conductivity at t=0 s, that is representative of the leader channel preceding the lightning stroke. It has been checked that the total energy in the initialization channel is very small compared to the total energy input during the lightning and that it has no significant influence on the results presented in this study.

For the modelling of composite materials and specifically the electrical conductivity, implementing in Taranis a tensorial electrical conductivity model was required. A modification in the equation 2 is done where the electrical conductivity is no longer a scalar variable $\sigma$ but a tensorial variable $\overline{\sigma}$. Its expression is given by equation 4.

$$ \overline{\overline{\sigma}} = \sigma_{\parallel} \cdot \overline{\overline{C}} = \sigma_{\parallel} \cdot \begin{bmatrix} \cos\theta^2 + \frac{\sigma_p}{\sigma_{\parallel}} \sin\theta^2 & \cos\theta \sin\theta \left[1 - \frac{\sigma_p}{\sigma_{\parallel}}\right] & 0 \\ \cos\theta \sin\theta \left[1 - \frac{\sigma_p}{\sigma_{\parallel}}\right] & \sin\theta^2 + \frac{\sigma_p}{\sigma_{\parallel}} \cos\theta^2 & 0 \\ 0 & 0 & \frac{\sigma_p}{\sigma_{\parallel}} \end{bmatrix} \quad (4) $$

The $\theta$ angle is formed between the x axis and the direction chosen for the composite fibers. In this work $\theta = 0$, and the x direction is chosen to be the direction of the carbon fibers and the y direction is chosen to be the direction perpendicular to the carbon fibers. The tensorial electrical conductivity depends on the value of the electrical conductivity along carbon fibers $\sigma_{\parallel}$, on the one perpendicular to the fibers $\sigma_{\perp}$ and the electrical conductivity in the direction transverse to the fibers through the material taken here in the z direction $\sigma_{\tan}$.
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

**3.2 Anisotropic mesh refinement workflow based on geflo.a**

Recent progress has matured anisotropic tetrahedra mesh adaptation to provide automatic control of estimated errors, which have been verified by code-to-code comparison. Alauzet documented the dramatic progress made during the last decade for anisotropic solution-adaptive methods to resolve simulations with shocks and boundary layers (Alauzet and Loseille, 2016). The adaptive computations use Hessian-based interpolation error estimates based on fully unstructured tetrahedra meshes. This method yields an anisotropic prescription of optimal meshes that are independent of geometric complexity. The mesh adaptation is based on the metric concept and the continuous mesh theory. The details of this framework are given in (Loseille and Alauzet, 2011a,b). Here, the main strategy is just recalled to tackle time-accurate mesh adaptation for unsteady simulations. An innovative mesh adaptation strategy for time-dependent problem based on a fixed-point algorithm has been developed in (Alauzet et al., 2018). The main goals are to guarantee a control of the spatial and temporal interpolation errors during the whole simulation and to reduce the number of mesh adaptations in order to decrease the error introduced by the transfers of solutions. This strategy starts with the observation that direct extension of steady adaptation algorithms to unsteady problems is not appropriate: specific algorithms must be developed that truly take into account the transient nature of the solution. It relies on the assumption that the temporal error is always controlled by the spatial one, which is indeed the case when solving, with an explicit time-integration, a linear advection problem under a CFL condition (Alauzet et al., 2007; Olivier, 2011). The fixed-point time-accurate anisotropic mesh adaptation algorithm relies on:

1. a *global* fixed-point mesh adaptation algorithm to converge the mesh adaptation non-linear problem,

2. the extension of the multi-scale error estimate (Loseille and Alauzet) to unsteady problems by proposing a $L^p$ space-time error analysis,

3. a conservative solution transfer operator to considerably diminish the error introduced by this stage of the algorithm (Alauzet, 2016),

4. a high-quality anisotropic local remeshing controlling the heights of the elements to guarantee optimal solver time-steps and to ensure maximal robustness in the meshing process. This issue has been handled with care in (Loseille, 2017).

The proposed approach overcomes all the problems relative to mesh adaptation for time-dependent simulations such as latency of the mesh with respect to the solution, time-accurate error estimate, and reduction of the error due to solution transfer from one mesh to another. In the present paper, the interpolation error is controlled by a scalar variable, namely $\rho E$ in the $L^2$ norm. In the course of the adaptation process, sequences of meshes that minimises the interpolation error in space and time are generated. This mesh adaptation framework is composed of a flow solver, here the MHD code Taranis, a metric-based error estimate, and a mesh adaptation code called geflo.a. The code geflo.a is being developed in Alauzet's team Gamma at INRIA France for more than ten years. A comprehensive and thorough description of the algorithm developed and used in geflo.a can be found in (Alauzet and Frazza, 2021) as well as in (Vanharen et al., 2022b) and (Vanharen et al., 2022a).

**3.3 Parameters of the 3D simulation and adaptation process**

Numerical simulation is carried out with a 3D configuration. The simulation domain is shown on Figure 1. The lightning arc is oriented in the z direction. The simulation domain size is $L_x = 40 \text{ cm} \times L_y = 40 \text{ cm} \times L_z = 30 \text{ cm}$. The thickness of the composite panel is $l_c=2 \text{ cm}$. The top domain above the composite panel is where the arc is simulated and is $l_f=15 \text{ cm}$ high. Boundary conditions for this simulated domain are detailed in Table 3.
It is here interesting to note here that a first mesh is generated by Gmsh ((Geuzaine and Remacle, 2009)) and consists of 291243 elements and 47744 nodes. For the workflow of the Taranis/geflo.a to take place, different parameters need to be valued: "$\rho E$" is the scalar variable followed for the error estimate in the feature base method, $h_{min} = 40 \ \mu\text{m}$ is the minimum length for an edge created in the mesh generated, $h_{max}= 0.5 \text{ m}$ is the maximum edge length desired and finally the limitation on the number of nodes generated in the new mesh that is equal here to 300 000 nodes. Computation was done on 96 MPI processes and took 10 days on the SATOR HPC machine at ONERA.
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

**3.2 Anisotropic mesh refinement workflow based on geflo.a**

Recent progress has matured anisotropic tetrahedra mesh adaptation to provide automatic control of estimated errors, which have been verified by code-to-code comparison. Alauzet documented the dramatic progress made during the last decade for anisotropic solution-adaptive methods to resolve simulations with shocks and boundary layers (Alauzet and Loseille, 2016). The adaptive computations use Hessian-based interpolation error estimates based on fully unstructured tetrahedra meshes. This method yields an anisotropic prescription of optimal meshes that are independent of geometric complexity. The mesh adaptation is based on the metric concept and the continuous mesh theory. The details of this framework are given in (Loseille and Alauzet, 2011a,b). Here, the main strategy is just recalled to tackle time-accurate mesh adaptation for unsteady simulations. An innovative mesh adaptation strategy for time-dependent problem based on a fixed-point algorithm has been developed in (Alauzet et al., 2018). The main goals are to guarantee a control of the spatial and temporal interpolation errors during the whole simulation and to reduce the number of mesh adaptations in order to decrease the error introduced by the transfers of solutions. This strategy starts with the observation that direct extension of steady adaptation algorithms to unsteady problems is not appropriate: specific algorithms must be developed that truly take into account the transient nature of the solution. It relies on the assumption that the temporal error is always controlled by the spatial one, which is indeed the case when solving, with an explicit time-integration, a linear advection problem under a CFL condition (Alauzet et al., 2007; Olivier, 2011). The fixed-point time-accurate anisotropic mesh adaptation algorithm relies on:

1. a *global* fixed-point mesh adaptation algorithm to converge the mesh adaptation non-linear problem,

2. the extension of the multi-scale error estimate (Loseille and Alauzet) to unsteady problems by proposing a $L^p$ space-time error analysis,

3. a conservative solution transfer operator to considerably diminish the error introduced by this stage of the algorithm (Alauzet, 2016),

4. a high-quality anisotropic local remeshing controlling the heights of the elements to guarantee optimal solver time-steps and to ensure maximal robustness in the meshing process. This issue has been handled with care in (Loseille, 2017).

The proposed approach overcomes all the problems relative to mesh adaptation for time-dependent simulations such as latency of the mesh with respect to the solution, time-accurate error estimate, and reduction of the error due to solution transfer from one mesh to another. In the present paper, the interpolation error is controlled by a scalar variable, namely $\rho E$ in the $L^2$ norm. In the course of the adaptation process, sequences of meshes that minimises the interpolation error in space and time are generated. This mesh adaptation framework is composed of a flow solver, here the MHD code Taranis, a metric-based error estimate, and a mesh adaptation code called geflo.a. The code geflo.a is being developed in Alauzet's team Gamma at INRIA France for more than ten years. A comprehensive and thorough description of the algorithm developed and used in geflo.a can be found in (Alauzet and Frazza, 2021) as well as in (Vanharen et al., 2022b) and (Vanharen et al., 2022a).

**3.3 Parameters of the 3D simulation and adaptation process**

Numerical simulation is carried out with a 3D configuration. The simulation domain is shown on Figure 1. The lightning arc is oriented in the z direction. The simulation domain size is $L_x = 40 \text{ cm} \times L_y = 40 \text{ cm} \times L_z = 30 \text{ cm}$. The thickness of the composite panel is $l_c=2 \text{ cm}$. The top domain above the composite panel is where the arc is simulated and is $l_f=15 \text{ cm}$ high. Boundary conditions for this simulated domain are detailed in Table 3.
It is here interesting to note here that a first mesh is generated by Gmsh ((Geuzaine and Remacle, 2009)) and consists of 291243 elements and 47744 nodes. For the workflow of the Taranis/geflo.a to take place, different parameters need to be valued: "$\rho E$" is the scalar variable followed for the error estimate in the feature base method, $h_{min} = 40 \ \mu\text{m}$ is the minimum length for an edge created in the mesh generated, $h_{max}= 0.5 \text{ m}$ is the maximum edge length desired and finally the limitation on the number of nodes generated in the new mesh that is equal here to 300 000 nodes. Computation was done on 96 MPI processes and took 10 days on the SATOR HPC machine at ONERA.
<figure><img src="image_2.png" /><figcaption>Figure 1: Simulation domain in 3D with numbered boundary surfaces.</figcaption></figure>

# 4. Simulation of an arc in a pulsed phase interacting with a composite panel

## 4.1 Experimental results of a pulsed arc interacting with a composite material

The goal of this study is to further the validation of the Taranis code that was done in 2D in (Pechereau et al., 2019). Here we compare the simulation results with one experiment that was carried out in the Ph.D. work of (Martins, 2016). Specifically in chapter 6 and page 130, it is shown the schematic, here on Figure 2, of an experiment where an arc in a pulsed phase is ignited above an unprotected composite material. By positioning two high speed cameras perpendicular to each other, it is observed that the root of the developing arc is forming in time a diamond shape. On Figure 3, it is interesting to observe that the arc root is developing faster in the direction perpendicular of the carbon fiber orientation compared to the parallel one. This behaviour can be explained by the fact that the electrical conductivity of the arc is higher compared to the composite in the perpendicular direction of propagation. In a previous work in (Pechereau et al., 2019), the electrical conductivity of an arc in pulsed phase driven by a "G" current wave was found roughly equal to $3.10^4$ S.m$^{-1}$. This value is very close to the one of carbon fibers equal to roughly $10^4$ S.m$^{-1}$. On the other hand in the direction perpendicular and tangential of the carbon fiber in the composite material, the electrical conductivity is lower by order of magnitudes. It was found in (Pechereau et al., 2019), that when the electrical conductivity of the arc is higher than the material on which the arc is forming, the current density is increasing at the edges of the arc root. By doing so, the Joule effect is also increasing and heats up the cold air around the arc. That effect coupled with radiative transfer that cools down the arc column and heats up also the edges of the arc are responsible for the spreading of the arc root close to the material surface. When the electrical conductivity of the material simulated is decreased from $10^4$ S.m$^{-1}$ to $10^2$ S.m$^{-1}$, the spreading of the arc root is increasing. In the experimental result shown on Figure 3, the arc root is spread on a wider surface in the direction perpendicular to the carbon fibers that has an electrical conductivity order of magnitude lower to the one parallel to the carbon fibers. Interestingly in the parallel direction of the carbon fibers, the arc root has not propagated further than the column above it. That means that the arc has an electrical conductivity in the range of the one of the carbon fibers. That experimental result is very interesting as it shows the influence of the tensorial electrical conductivity of a composite material on the dynamics of an arc in a pulsed phase but also provide experimental results that can be compared with simulation ones.
<figure><img src="image_3.png" /><figcaption>Figure 1: Simulation domain in 3D with numbered boundary surfaces.</figcaption></figure>

# 4. Simulation of an arc in a pulsed phase interacting with a composite panel

## 4.1 Experimental results of a pulsed arc interacting with a composite material

The goal of this study is to further the validation of the Taranis code that was done in 2D in (Pechereau et al., 2019). Here we compare the simulation results with one experiment that was carried out in the Ph.D. work of (Martins, 2016). Specifically in chapter 6 and page 130, it is shown the schematic, here on Figure 2, of an experiment where an arc in a pulsed phase is ignited above an unprotected composite material. By positioning two high speed cameras perpendicular to each other, it is observed that the root of the developing arc is forming in time a diamond shape. On Figure 3, it is interesting to observe that the arc root is developing faster in the direction perpendicular of the carbon fiber orientation compared to the parallel one. This behaviour can be explained by the fact that the electrical conductivity of the arc is higher compared to the composite in the perpendicular direction of propagation. In a previous work in (Pechereau et al., 2019), the electrical conductivity of an arc in pulsed phase driven by a "G" current wave was found roughly equal to $3.10^4$ S.m$^{-1}$. This value is very close to the one of carbon fibers equal to roughly $10^4$ S.m$^{-1}$. On the other hand in the direction perpendicular and tangential of the carbon fiber in the composite material, the electrical conductivity is lower by order of magnitudes. It was found in (Pechereau et al., 2019), that when the electrical conductivity of the arc is higher than the material on which the arc is forming, the current density is increasing at the edges of the arc root. By doing so, the Joule effect is also increasing and heats up the cold air around the arc. That effect coupled with radiative transfer that cools down the arc column and heats up also the edges of the arc are responsible for the spreading of the arc root close to the material surface. When the electrical conductivity of the material simulated is decreased from $10^4$ S.m$^{-1}$ to $10^2$ S.m$^{-1}$, the spreading of the arc root is increasing. In the experimental result shown on Figure 3, the arc root is spread on a wider surface in the direction perpendicular to the carbon fibers that has an electrical conductivity order of magnitude lower to the one parallel to the carbon fibers. Interestingly in the parallel direction of the carbon fibers, the arc root has not propagated further than the column above it. That means that the arc has an electrical conductivity in the range of the one of the carbon fibers. That experimental result is very interesting as it shows the influence of the tensorial electrical conductivity of a composite material on the dynamics of an arc in a pulsed phase but also provide experimental results that can be compared with simulation ones.
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

<table><caption>Table 3: Boundary conditions for the arc simulation. N: Neumann BC; ø: input voltage on BC n°2; S.O: Supersonic Outflow; A.S.W: Adiabatic Slip Wall.</caption><thead><tr><th>Variable</th><th>BC n°1</th><th>BC n°2</th><th>BC n°3</th><th>BC n°4</th><th>BC n°5</th><th>BC n°6</th><th>BC n°7</th><th>BC n°8</th></tr></thead><tbody><tr><td>φ</td><td>0</td><td>ø</td><td>0</td><td>N</td><td>0</td><td>0</td><td>N</td><td>0</td></tr><tr><td>Euler</td><td>/</td><td>A.S.W</td><td>/</td><td>S.O</td><td>/</td><td>/</td><td>S.O</td><td>/</td></tr><tr><td>A<sub>x</sub></td><td>0</td><td>0</td><td>N</td><td>N</td><td>0</td><td>N</td><td>N</td><td>0</td></tr><tr><td>A<sub>y</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>N</td><td>0</td><td>0</td><td>N</td></tr><tr><td>A<sub>z</sub></td><td>N</td><td>N</td><td>0</td><td>0</td><td>N</td><td>0</td><td>0</td><td>0</td></tr><tr><td>G<sub>v</sub></td><td>0</td><td>N</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table><table><thead><tr><th>Variable</th><th>BC n°9</th><th>BC n°10</th><th>BC n°11</th><th>BC n°12</th><th>BC n°13</th><th>BC n°14</th><th>BC n°15</th><th>BC n°16</th></tr></thead><tbody><tr><td>φ</td><td>N</td><td>N</td><td>N</td><td>N</td><td>N</td><td>/</td><td>/</td><td>N</td></tr><tr><td>Euler</td><td>S.O</td><td>/</td><td>/</td><td>/</td><td>/</td><td>A.S.W</td><td>/</td><td>S.O</td></tr><tr><td>A<sub>x</sub></td><td>0</td><td>0</td><td>N</td><td>0</td><td>N</td><td>/</td><td>/</td><td>0</td></tr><tr><td>A<sub>y</sub></td><td>N</td><td>N</td><td>0</td><td>N</td><td>0</td><td>/</td><td>/</td><td>N</td></tr><tr><td>A<sub>z</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>/</td><td>/</td><td>0</td></tr><tr><td>G<sub>v</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>/</td><td>/</td><td>0</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic diagram of the different camera positions. Chapter 6, page 130 in (Martins, 2016)</figcaption></figure>

## 4.2 Simulation results in 3D using anisotropic mesh refinement

In order to simulate the effect of conductivity of the composite panel on the arc dynamics, a tensorial conductivity is used to model the conductivity in the composite panel and requires to do 3D computation. For the purpose of the qualitative comparison with experiments, order of magnitude values were used as a first step for the electrical conductivity. Following the work of (Jennings and Hardwick, 1992) and (Lago, 2004), we used these values for the tensorial conductivity of the single panel simulated in this work: $\sigma_{\parallel}=10^4$ S.m$^{-1}$, $\sigma_{\perp}=10^2$ S.m$^{-1}$ and $\sigma_{tan}=10^2$ S.m$^{-1}$. The orientation of the carbon fibers in the composite is set in the x direction and so is the direction taken for setting $\sigma_{\parallel}$. In order to follow in 3D the dynamics of the arc under a "G" current wave, an isocontour of temperature T = 5 000 K is presented on Figure 4 at different times t=2 µs, 10 µs and 20 µs. As in the experiment, the arc root is spreading farther in the direction perpendicular to the carbon fibers along the y direction compared to its spread in the direction of the carbon fibers set in the x direction. Starting from a uniform column, in time the spreading of the arc root is clearly forming a diamond shape as in the experiment. In this simulation, aside from these 3D results, a workflow that couples the MHD code Taranis and the anisotropic mesh refinement code geflo.a was used. Specifically the unsteady workflow was used. Then during this simulation the mesh followed in time the development of the arc simulated. Figure 5 shows at t=2 µs, 2D cuts along the y direction on top of the image and along the x direction on the bottom of the image. On the left part of the image 2D distribution of temperature is displayed with the corresponding mesh on the right. We observe that early in the simulation, although we observe already a different spreading of the arc root in the two perpendicular directions, the anisotropic mesh generated along the two directions is quite similar. On the other hand,
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

<table><caption>Table 3: Boundary conditions for the arc simulation. N: Neumann BC; ø: input voltage on BC n°2; S.O: Supersonic Outflow; A.S.W: Adiabatic Slip Wall.</caption><thead><tr><th>Variable</th><th>BC n°1</th><th>BC n°2</th><th>BC n°3</th><th>BC n°4</th><th>BC n°5</th><th>BC n°6</th><th>BC n°7</th><th>BC n°8</th></tr></thead><tbody><tr><td>φ</td><td>0</td><td>ø</td><td>0</td><td>N</td><td>0</td><td>0</td><td>N</td><td>0</td></tr><tr><td>Euler</td><td>/</td><td>A.S.W</td><td>/</td><td>S.O</td><td>/</td><td>/</td><td>S.O</td><td>/</td></tr><tr><td>A<sub>x</sub></td><td>0</td><td>0</td><td>N</td><td>N</td><td>0</td><td>N</td><td>N</td><td>0</td></tr><tr><td>A<sub>y</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>N</td><td>0</td><td>0</td><td>N</td></tr><tr><td>A<sub>z</sub></td><td>N</td><td>N</td><td>0</td><td>0</td><td>N</td><td>0</td><td>0</td><td>0</td></tr><tr><td>G<sub>v</sub></td><td>0</td><td>N</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

<table><thead><tr><th>Variable</th><th>BC n°9</th><th>BC n°10</th><th>BC n°11</th><th>BC n°12</th><th>BC n°13</th><th>BC n°14</th><th>BC n°15</th><th>BC n°16</th></tr></thead><tbody><tr><td>φ</td><td>N</td><td>N</td><td>N</td><td>N</td><td>N</td><td>/</td><td>/</td><td>N</td></tr><tr><td>Euler</td><td>S.O</td><td>/</td><td>/</td><td>/</td><td>/</td><td>A.S.W</td><td>/</td><td>S.O</td></tr><tr><td>A<sub>x</sub></td><td>0</td><td>0</td><td>N</td><td>0</td><td>N</td><td>/</td><td>/</td><td>0</td></tr><tr><td>A<sub>y</sub></td><td>N</td><td>N</td><td>0</td><td>N</td><td>0</td><td>/</td><td>/</td><td>N</td></tr><tr><td>A<sub>z</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>/</td><td>/</td><td>0</td></tr><tr><td>G<sub>v</sub></td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>/</td><td>/</td><td>0</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 2: Schematic diagram of the different camera positions. Chapter 6, page 130 in (Martins, 2016)</figcaption></figure>

4.2 Simulation results in 3D using anisotropic mesh refinement

In order to simulate the effect of conductivity of the composite panel on the arc dynamics, a tensorial conductivity is used to model the conductivity in the composite panel and requires to do 3D computation. For the purpose of the qualitative comparison with experiments, order of magnitude values were used as a first step for the electrical conductivity. Following the work of (Jennings and Hardwick, 1992) and (Lago, 2004), we used these values for the tensorial conductivity of the single panel simulated in this work: $\sigma_{||}=10^4$ S.m$^{-1}$, $\sigma_{\perp}=10^2$ S.m$^{-1}$ and $\sigma_{tan}=10^2$ S.m$^{-1}$. The orientation of the carbon fibers in the composite is set in the x direction and so is the direction taken for setting $\sigma_{||}$. In order to follow in 3D the dynamics of the arc under a "G" current wave, an isocontour of temperature T = 5 000 K is presented on Figure 4 at different times t=2 µs, 10 µs and 20 µs. As in the experiment, the arc root is spreading farther in the direction perpendicular to the carbon fibers along the y direction compared to its spread in the direction of the carbon fibers set in the x direction. Starting from a uniform column, in time the spreading of the arc root is clearly forming a diamond shape as in the experiment. In this simulation, aside from these 3D results, a workflow that couples the MHD code Taranis and the anisotropic mesh refinement code geflo.a was used. Specifically the unsteady workflow was used. Then during this simulation the mesh followed in time the development of the arc simulated. Figure 5 shows at t=2 µs, 2D cuts along the y direction on top of the image and along the x direction on the bottom of the image. On the left part of the image 2D distribution of temperature is displayed with the corresponding mesh on the right. We observe that early in the simulation, although we observe already a different spreading of the arc root in the two perpendicular directions, the anisotropic mesh generated along the two directions is quite similar. On the other hand,
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

<figure><img src="image_4.png" /><figcaption>Figure 3: Pictures of the arc root evolution over time applied to a carbon fiber composite panel. (a) Camera axis perpendicular to the carbon fiber orientation. (b) Camera axis parallel to the carbon fiber orientation. Taken from (Martins, 2016)</figcaption></figure>

at 10 µs and 20 µs, on Figure 6 and 7, we clearly see the difference in the mesh used along the two directions and that mesh points are placed precisely at the locations that outlines the shape and the root of the spreading arc.

It is important to point out that without mesh refinement, that simulation in 3D would not have been possible with the resources used for this work. With a fixed mesh and a mesh size of 40 µm, the number of cells required is in the order of $10^9-10^{10}$ cells. That simulation would have required more than 10 000 MPI processes for at least a few weeks without taking into account the scalability of the Taranis code and specifically the scalability of the linear solvers used in it. Compared to the simulation presented here that was carried out on 96 MPI processes and using no more than 300 000 nodes or around $1.8 \ 10^6$ cells, the order of magnitude of savings in computing power is greater than 2.

The results presented in this work are twofold: the 3D simulation using the workflow Taranis/feflo.a is successful in recreating the diamond shape of the arc root on a composite material and is also the first time that anisotropic mesh refinement is used with a MHD code to simulate the dynamics of an arc in a pulsed phase in interaction with a composite material.

<figure><img src="image_5.png" /><figcaption>Figure 4: 3D plots of isocontours value at T = 5 000 K at times t=2 µs, 10 µs and 20 µs are shown.</figcaption></figure>
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

<figure><img src="image_5.png" /><figcaption>Figure 3: Pictures of the arc root evolution over time applied to a carbon fiber composite panel. (a) Camera axis perpendicular to the carbon fiber orientation. (b) Camera axis parallel to the carbon fiber orientation. Taken from (Martins, 2016)</figcaption></figure>

at 10 µs and 20 µs, on Figure 6 and 7, we clearly see the difference in the mesh used along the two directions and that mesh points are placed precisely at the locations that outlines the shape and the root of the spreading arc.

It is important to point out that without mesh refinement, that simulation in 3D would not have been possible with the resources used for this work. With a fixed mesh and a mesh size of 40 µm, the number of cells required is in the order of $10^9-10^{10}$ cells. That simulation would have required more than 10 000 MPI processes for at least a few weeks without taking into account the scalability of the Taranis code and specifically the scalability of the linear solvers used in it. Compared to the simulation presented here that was carried out on 96 MPI processes and using no more than 300 000 nodes or around $1.8 \ 10^6$ cells, the order of magnitude of savings in computing power is greater than 2.

The results presented in this work are twofold: the 3D simulation using the workflow Taranis/eflo.a is successful in recreating the diamond shape of the arc root on a composite material and is also the first time that anisotropic mesh refinement is used with a MHD code to simulate the dynamics of an arc in a pulsed phase in interaction with a composite material.

<figure><img src="image_6.png" /><figcaption>Figure 4: 3D plots of isocontours value at T = 5 000 K at times t=2 µs, 10 µs and 20 µs are shown.</figcaption></figure>
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

<figure><img src="image_6.png" /><figcaption>Figure 5: 2D plots of temperature distributions and mesh at times t = 2 µs. 2D cuts along the y direction "(a)" and along the x direction "(b)"</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: 2D plots of temperature distributions and mesh at times t = 10 µs. 2D cuts along the y direction "(a)" and along the x direction "(b)"</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: 2D plots of temperature distributions and mesh at times t = 20 µs. 2D cuts along the y direction "(a)" and along the x direction "(b)"</figcaption></figure>

**5. Conclusions**

A new workflow has been developed and used to study into details the dynamics of an arc in a pulsed phase in in-
teraction with a composite material. Thanks to the coupling of the 3D MHD code Taranis with the anisotropic mesh
refinement code geflo.a, savings in computing power is greater than 2 order of magnitudes and allowed this simulation
to be done on 96 MPI processes and 10 days of computations. In doing so, a qualitative comparison was achieved and
the simulation showed, as in the experiments, the creation of a diamond shape by the arc root of an arc in a pulsed
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

<figure><img src="image_7.png" /><figcaption>Figure 5: 2D plots of temperature distributions and mesh at times t = 2 µs. 2D cuts along the y direction "(a)" and along the x direction "(b)"</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: 2D plots of temperature distributions and mesh at times t = 10 µs. 2D cuts along the y direction "(a)" and along the x direction "(b)"</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: 2D plots of temperature distributions and mesh at times t = 20 µs. 2D cuts along the y direction "(a)" and along the x direction "(b)"</figcaption></figure>

**5. Conclusions**

A new workflow has been developed and used to study into details the dynamics of an arc in a pulsed phase in in-
teraction with a composite material. Thanks to the coupling of the 3D MHD code Taranis with the anisotropic mesh
refinement code geflo.a, savings in computing power is greater than 2 order of magnitudes and allowed this simulation
to be done on 96 MPI processes and 10 days of computations. In doing so, a qualitative comparison was achieved and
the simulation showed, as in the experiments, the creation of a diamond shape by the arc root of an arc in a pulsed
phase propagating on an unprotected composite material.

Based on these results, more work will be dedicated on modelling more realistically the composite material used in the experiment. The first step is to decrease the composite panel thickness from to 2 cm to 2mm as in the experiments. The second step is to take into account the 8 layers that composes the composite material with each layer having different values of tensorial electrical conductivity.

# 6. Acknowledgments

The authors wish to thank the French Civil Aviation Authority (DGAC), France Relance and NextGenerationEU for their supports. The authors would like to also thank Onera for the use of the supercomputer SATOR.

# References

F. Alauzet. A parallel matrix-free conservative solution interpolation on unstructured tetrahedral meshes. *Comput. Methods Appl. Mech. Engrg.*, 299:116–142, 2016.

F. Alauzet and L. Frazza. Feature-based and goal-oriented anisotropic mesh adaptation for rans applications in aero-nautics and aerospace. *Journal of Computational Physics*, 439:110340, 2021.

F. Alauzet and A. Loseille. A decade of progress on anisotropic mesh adaptation for computational fluid dynamics. *Computer-Aided Design*, 72:13–39, 2016.

F. Alauzet, P. Frey, P. George, and B. Mohammadi. 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations. *J. Comp. Phys.*, 222:592–623, 2007.

F. Alauzet, A. Loseille, and G. Olivier. Time-accurate multi-scale anisotropic mesh adaptation for unsteady flows in CFD. *J. Comp. Phys.*, 373:28–63, 2018.

S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith. Efficient management of parallelism in object oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, *Modern Software Tools in Scientific Computing*, pages 163–202. Birkhäuser Press, 1997.

S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang. PETSc/TAO users manual. Technical Report ANL-21/39 - Revision 3.19, Argonne National Laboratory, 2023a.

S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang. PETSc Web page. https://petsc.org/, 2023b. URL https://petsc.org/.

S. Chauveau, C. Deron, M.-Y. Perrin, P. Rivière, and A. Soufiani. Radiative transfer in LTE air plasmas for temperatures up to 15,000 k. *Journal of Quantitative Spectroscopy and Radiative Transfer*, 77(2):113–130, 2003.

A. D'angola, G. Colonna, C. Gorse, and M. Capitelli. Thermodynamic and transport properties in equilibrium air plasmas in a wide pressure and temperature range. *The European Physical Journal D*, 46(1):129–150, 2008.

F. A. Fisher and J. A. Plumer. Lightning protection of aircraft. Technical report, 1977.

C. Geuzaine and J.-F. Remacle. Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities. *International journal for numerical methods in engineering*, 79(11):1309–1331, 2009.

N. Jennings and C. Hardwick. A computational approach to predicting the extent of arc root damage in cfc panels, 15th int. in *Aerospace and Ground Conf. on Lightning and Static Electricity (Atlantic City) pp*, pages 41–1, 1992.

F. Lago. Modélisation de l'interaction entre un arc électrique et une surface: application au foudroyement d'un aéronef. PhD thesis, Toulouse 3, 2004.
phase propagating on an unprotected composite material.

Based on these results, more work will be dedicated on modelling more realistically the composite material used in the experiment. The first step is to decrease the composite panel thickness from to 2 cm to 2mm as in the experiments. The second step is to take into account the 8 layers that composes the composite material with each layer having different values of tensorial electrical conductivity.

# 6. Acknowledgments

The authors wish to thank the French Civil Aviation Authority (DGAC), France Relance and NextGenerationEU for their supports. The authors would like to also thank Onera for the use of the supercomputer SATOR.

# References

F. Alauzet. A parallel matrix-free conservative solution interpolation on unstructured tetrahedral meshes. *Comput. Methods Appl. Mech. Engrg.*, 299:116–142, 2016.

F. Alauzet and L. Frazza. Feature-based and goal-oriented anisotropic mesh adaptation for rans applications in aero-nautics and aerospace. *Journal of Computational Physics*, 439:110340, 2021.

F. Alauzet and A. Loseille. A decade of progress on anisotropic mesh adaptation for computational fluid dynamics. *Computer-Aided Design*, 72:13–39, 2016.

F. Alauzet, P. Frey, P. George, and B. Mohammadi. 3D transient fixed point mesh adaptation for time-dependent problems: Application to CFD simulations. *J. Comp. Phys.*, 222:592–623, 2007.

F. Alauzet, A. Loseille, and G. Olivier. Time-accurate multi-scale anisotropic mesh adaptation for unsteady flows in CFD. *J. Comp. Phys.*, 373:28–63, 2018.

S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith. Efficient management of parallelism in object oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, *Modern Software Tools in Scientific Computing*, pages 163–202. Birkhäuser Press, 1997.

S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang. PETSc/TAO users manual. Technical Report ANL-21/39 - Revision 3.19, Argonne National Laboratory, 2023a.

S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang. PETSc Web page. https://petsc.org/, 2023b. URL https://petsc.org/.

S. Chauveau, C. Deron, M.-Y. Perrin, P. Rivière, and A. Soufiani. Radiative transfer in lte air plasmas for temperatures up to 15,000 k. *Journal of Quantitative Spectroscopy and Radiative Transfer*, 77(2):113–130, 2003.

A. D'angola, G. Colonna, C. Gorse, and M. Capitelli. Thermodynamic and transport properties in equilibrium air plasmas in a wide pressure and temperature range. *The European Physical Journal D*, 46(1):129–150, 2008.

F. A. Fisher and J. A. Plumer. Lightning protection of aircraft. Technical report, 1977.

C. Geuzaine and J.-F. Remacle. Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities. *International journal for numerical methods in engineering*, 79(11):1309–1331, 2009.

N. Jennings and C. Hardwick. A computational approach to predicting the extent of arc root damage in cfc panels, 15th int. in *Aerospace and Ground Conf. on Lightning and Static Electricity (Atlantic City) pp*, pages 41–1, 1992.

F. Lago. *Modélisation de l'interaction entre un arc électrique et une surface: application au foudroyement d'un aéronef*. PhD thesis, Toulouse 3, 2004.
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

A. Loseille. Chapter 10 - Unstructured Mesh Generation and Adaptation. In R. Abgrall and C.-W. Shu, editors, *Handbook of Numerical Methods for Hyperbolic Problems*, volume 18 of *Handbook of Numerical Analysis*, pages 263–302. Elsevier, 2017.

A. Loseille and F. Alauzet. Optimal 3D highly anisotropic mesh adaptation based on the continuous mesh framework. In *P*.

A. Loseille and F. Alauzet. Continuous mesh framework. Part I: well-posed continuous interpolation error. *SIAM J. Numer. Anal.*, 49(1):38–60, 2011a.

A. Loseille and F. Alauzet. Continuous mesh framework. Part II: validations and applications. *SIAM J. Numer. Anal.*, 49(1):61–86, 2011b.

R. A. S. Martins. *Etude expérimentale et théorique d'un arc de foudre et son interaction avec un matériau aéronautique*. PhD thesis, Université Paris-Saclay, 2016.

G. Olivier. *Anisotropic metric-based mesh adaptation for unsteady CFD simulations involving moving geometries*. PhD thesis, Université Pierre et Marie Curie, Paris VI, Paris, France, 2011.

F. Pechereau, F. Tholin, C. Zaepffel, R. Sousa Martins, and P. Lalande. Mhd simulation of a lightning column during the pulsed arc phase. 2019.

F. Tholin, L. Chemartin, and P. Lalande. Numerical investigation of the interaction of a lightning and an aeronautic skin during the pulsed arc phase. 2015.

F. Tholin, L. Chemartin, and P. Lalande. Simulation of the lightning arc root interaction with anisotropic materials. In *ICOLSE 2017*, 2017.

J. Vanharen, A. Loseille, and F. Alauzet. Non-manifold anisotropic mesh adaptation: application to fluid-structure interaction. *Engineering with Computers*, 38(5):4269–4288, 2022a.

J. Vanharen, A. Loseille, F. Alauzet, and M. A. Park. Near-field anisotropic mesh adaptation for the third AIAA sonic boom workshop. *Journal of Aircraft*, 59(3):683–696, 2022b.
3D MODELING OF AN ARC IN INTERACTION WITH A COMPOSITE PANEL

A. Loseille. Chapter 10 - Unstructured Mesh Generation and Adaptation. In R. Abgrall and C.-W. Shu, editors, *Handbook of Numerical Methods for Hyperbolic Problems*, volume 18 of *Handbook of Numerical Analysis*, pages 263–302. Elsevier, 2017.

A. Loseille and F. Alauzet. Optimal 3D highly anisotropic mesh adaptation based on the continuous mesh framework. In *P*.

A. Loseille and F. Alauzet. Continuous mesh framework. Part I: well-posed continuous interpolation error. *SIAM J. Numer. Anal.*, 49(1):38–60, 2011a.

A. Loseille and F. Alauzet. Continuous mesh framework. Part II: validations and applications. *SIAM J. Numer. Anal.*, 49(1):61–86, 2011b.

R. A. S. Martins. *Etude expérimentale et théorique d'un arc de foudre et son interaction avec un matériau aéronautique*. PhD thesis, Université Paris-Saclay, 2016.

G. Olivier. *Anisotropic metric-based mesh adaptation for unsteady CFD simulations involving moving geometries*. PhD thesis, Université Pierre et Marie Curie, Paris VI, Paris, France, 2011.

F. Pechereau, F. Tholin, C. Zaepffel, R. Sousa Martins, and P. Lalande. Mhd simulation of a lightning column during the pulsed arc phase. 2019.

F. Tholin, L. Chemartin, and P. Lalande. Numerical investigation of the interaction of a lightning and an aeronautic skin during the pulsed arc phase. 2015.

F. Tholin, L. Chemartin, and P. Lalande. Simulation of the lightning arc root interaction with anisotropic materials. In *ICOLSE 2017*, 2017.

J. Vanharen, A. Loseille, and F. Alauzet. Non-manifold anisotropic mesh adaptation: application to fluid-structure interaction. *Engineering with Computers*, 38(5):4269–4288, 2022a.

J. Vanharen, A. Loseille, F. Alauzet, and M. A. Park. Near-field anisotropic mesh adaptation for the third AIAA sonic boom workshop. *Journal of Aircraft*, 59(3):683–696, 2022b.