Fluid–rigid body interaction in a compressible electrically conducting fluid

We consider a system of multiple insulating rigid bodies moving inside of an electrically conducting compressible fluid. In this system, we take into account the interaction of the fluid with the bodies as well as with the electromagnetic fields trespassing both the fluid and the solids. The main result of this paper yields the existence of weak solutions to the system. While the mechanical part of the problem can be dealt with via a classical penalization method, the electromagnetic part requires an approximation by means of a hybrid discrete–continuous in time system: The discrete part of the approximation enables us to handle the solution‐dependent test functions in our variational formulation of the induction equation, whereas the continuous part makes sure that the nonnegativity of the density and subsequently a meaningful energy inequality is preserved in the approximate system.


Introduction
The goal of this article is the proof of the existence of weak solutions to a system of PDEs modelling the motion of multiple non-conducting rigid objects inside of an electrically conducting compressible fluid taking the interplay with the electromagnetic fields in these materials into account.The problem can be regarded as belonging to both the research areas of magnetohydrodynamics (MHD) and fluidstructure interaction (FSI).Indeed, MHD (c.f.[5,8,35]), on the one hand, describes a coupling between the Maxwell system and the Navier-Stokes equations which models the influence of the electromagnetic fields on the electrically conducting fluid and vice versa.FSI, on the other hand, models the interaction between the fluid and the solid bodies through a coupling between the Navier-Stokes system and the balance of mass and momentum of the rigid bodies.The interplay between the electromagnetic fields and the solids occurs implicitly via the fluid due to the non-conductivity of the solid material.This paper represents the extension of [1], where the corresponding problem was studied and solved for incompressible fluids, to the compressible case.Potential applications can be found in the area of biomechanics.A specific example is constituted by capsule endoscopy, a medical procedure during which small cameradevices are sent through the body with the aim of detecting diseases.Such endocapsules of a microscopic scale can be navigated through the (electrically conducting) blood by controlling their movement via the application of electromagnetic forces, c.f. [6,Section 4.4], [4].This technique can also be applied in the problem of drug delivery, in which microrobots are constructed to deliver drugs directly to the targeted area of the body without affecting healthy tissue.However, in view of the usage of electrically conducting microrobots for these purposes an extension of the current model to electrically conducting rigid bodies would be needed.In order to classify our result within the wide range of works in MHD and FSI, we give a brief summary of the related literature.The existence of weak solutions to the MHD problem for a compressible fluid but without any rigid bodies involved is for example shown in [41] and [2].The fluids considered in the latter one of these articles are electrically as well as thermally conducting and moreover the existence of strong solutions is addressed therein.The existence of weak solutions to the MHD problem in the incompressible case is proved in [25].On the FSI side of the problem we first mention [23,42] for a general introduction to the fluid-rigid body interaction problem.Early existence results for weak solutions to this problem were obtained mainly in the incompressible case and include the articles [7,9,29,32], wherein the proof of local in time existence, i.e. existence up to a contact between several bodies or a body and the domain boundary, is achieved in the case of two and three spatial dimensions.The corresponding problem in the compressible case was considered in [10].A proof of the global in time existence of weak solutions to a model describing the interaction between multiple rigid bodies and a compressible fluid is given in [17].Corresponding results are also available for incompressible fluids, c.f. [40] for the twodimensional and [18] for the three-dimensional case.The article [40] in particular touches upon the question whether contacts between the bodies with each other or the domain boundary are possible and shows that such collisions can only occur if the relative acceleration and velocity between the colliding objects vanish.Moreover, the problem has been studied in the case of the Navier-slip instead of the classical no-slip boundary condition, c.f. [36], and also the question about the existence of strong solutions has been investigated in both the compressible, c.f. [3,30,31,39], and the incompressible case, c.f. [24,44,45].The articles [27] (for the case of two spatial dimensions) and [28] (for the 3D case) can be considered as a first step towards the coupling between the MHD and the FSI problem.The authors thereof studied a model of an incompressible electrically conducting fluid flowing around a fixed non-conducting solid region.This model served as the basis for the electromagnetic part of the model in [1], wherein Benešová, Nečasová, Schlömerkemper and the author of this article showed the (local in time) existence of weak solutions to the problem of one movable insulating rigid object travelling through an electrically conducting incompressible fluid.The main goal of the present article is to extend the latter result to the case of a compressible fluid.More precisely, we are able to prove the global in time existence of weak solutions to the interaction problem of finitely many insulating rigid bodies, an electrically conducting compressible fluid and the electromagnetic fields trespassing these materials, see Section 2.3.As in the incompressible case in [1], the main difficulty in the proof of the existence result is caused by the test functions in the weak formulation of the induction equation, c.f. ( 18), (29) below, which are chosen such that they depend on the solid region.While such test functions do not generate any problems in the case of an immovable solid region (see e.g. the proofs of [27, Theorem 2.1] and [28,Theorem 2.3]), difficulties arise in our scenario, where the solid domain depends on the overall solution to the system, which causes our problem to be highly coupled.Following the proof in the incompressible situation, we thus make use of a time discretization, which allows us to deal with this problem by decoupling it: At each discrete time we first calculate the domain of the solid bodies, which suggests a suitable definition for the test functions in the induction equation at that specific time.Only after this we solve the induction equation itself, which can then be achieved via standard methods.In the compressible situation, however, this procedure turns out to cause more problems.In particular, the author could not find a suitable way to discretize the Navier-Stokes system while preserving the non-negativity of the density.This is essential to obtain the uniform bounds from the energy inequality required to pass to the limit in the approximate system.We handle this problem by choosing a hybrid approximation system, in which the induction equation is discretized in time via the Rothe method ([38, Section 8]), whereas the mechanical part of the system is studied as a time-dependent problem on the small intervals between the discrete time points.The non-negativity of the density can then be derived by classical arguments and, by choosing the coupling terms in a suitable way, the discrete electromagnetic part and the continuous mechanical part of the system can be combined to an energy inequality with all desired features.To a smaller extent, we already used such a hybrid approximation in the incompressible situation [1], where, however, the time dependency was restricted to the transport equation for the characteristic function of the solid region.The expansion of this idea to the whole mechanical part of the system, in order to deal with the problems outlined above, is what constitutes the main novelty in the proof of our main result.The problem of the solution dependent test functions also appears in the weak formulation of the momentum equation, c.f. ( 17), (28) below.In this situation, however, we have the penalization method used e.g. in [17] and [40] readily available which allows us to evade the problem.More precisely, in this penalization method a sequence of approximate solutions to some fluid-only problems with classical test functions is constructed.Passing to the limit in this approximation one then returns to a fluidrigid body interaction system by letting the viscosity of the fluid rise to infinity in the later solid region.The paper is structured as follows: We begin by summarizing the notation needed for our model in Section 1.1 and subsequently present the model, divided into a mechanical and an electromagnetic part, in Sections 1.2 and 1.3.After introducing some additional notation in Section 2.1, we present the variational formulation of the above model in Section 2.2 as well as our main result in Section 2.3.In Section 3 we explain the main ideas for the proof of this result, which is based on an approximation of the original system.In Section 4 we solve this approximate problem and finally, in the remaining Sections 5-9, we pass to the limit in the approximation, proving the existence of a weak solution to the original system.

Model
The model we consider describes several non-conducting rigid bodies travelling through an electrically conducting compressible fluid as well as the involved electromagnetic fields.This model is a combination of (i) the mechanical fluid-rigid body interaction model used in [17] and (ii) the Maxwell system in the model used in [27,28] for the description of an electrically conducting fluid surrounding an immovable, non-conducting solid region.It is further the extension of the corresponding model for the incompressible case in [1] to the compressible situation.Let T ą 0 and let Ω Ă R 3 be a bounded domain.Inside of Ω we consider m P N insulating rigid bodies, the positions of which at time t P r0, T s are described through subsets S i ptq Ă Ω, i " 1, ..., m.The complement of the solid domain, F ptq :" ΩzSptq, with Sptq :" contains an electrically conducting viscous nonhomogeneous compressible fluid.We denote by Q the time-space domain Q :" p0, T q ˆΩ, which we split into a solid part Q s and a fluid part For any function defined on Q we mark its restriction to Q f or Q s by the superscript f or s, respectively.The interaction between the fluid, the solids and the electromagnetic fields in the domain is characterized through the mass density ρ : Q Ñ R, the velocity field u : Q Ñ R 3 , the magnetic induction B : Q Ñ R 3 , the electric field E : Q Ñ R 3 and the electric current j : Q Ñ R 3 .As indicated above, our overall model, which determines these functions, can be divided into two subsystems: The mechanical subsystem for the description of ρ and u and the electromagnetic subsystem for the description of B, E and j.

Mechanical subsystem
The mechanical quantities evolve according to the compressible Navier-Stokes equations in the fluid domain and the balance of linear and angular momentum of the rigid bodies in the solid region, respectively, rT ´pIds ¨n dσ `żS i ptq ρg dx, t P r0, T s, i " 1, ..., m, (4) px ´hi q ˆrT ´pIds n dσ `żS i ptq ρ `x ´hi ˘ˆg dx, t P r0, T s, i " 1, ..., m, (5) combined with the boundary and interface conditions uptq " 0 on BΩ, u f ptq ´us ptq " 0 on BSptq.
The identity ( 2) is known as the continuity equation.In the momentum equation ( 3) we see the pressure p, for which we assume the isentropic constitutive relation and the stress tensor T " Tpuq :" 2νDpuq `λId div u, Dpuq :" with the viscosity coefficients ν, λ P R which satisfy ν ą 0, λ `ν ě 0.Moreover, we have two forcing terms: The external force g and the Lorentz force 1 µ curlB ˆB with the magnetic permeability µ ą 0. We assume that µ is constant in the whole domain Q.
Since, in general, the magnetic permeability takes different values in conducting and insulating materials, this assumption is physically restrictive.However, it is necessary for the transition conditions on the magnetic induction B, c.f. Section 1.3 below.The Lorentz force constitutes the connection of the mechanical to the electromagnetic subsystem.In the balance of mass (4) and momentum (5) of the rigid bodies instead the Lorentz force does not appear, which is in accordance with the assumption of the bodies being non-conducting.The notation used in these relations includes the total mass m i of the i-th body, its center of mass h i and the associated inertia tensor J i , ρpt, xqx dx, t P r0, T s, The equations ( 4) and (5) determine the translational and rotational velocities V i and w i of the i-th rigid body, respectively, allowing us to express its total velocity as a rigid velocity field upt, xq " u s i pt, xq :" V i ptq `wi ptq ˆ`x ´hi ptq ˘for t P r0, T s, x P S i ptq.
The coupling between the fluid and the solids is incorporated into the surface integrals in ( 4) and ( 5) and into the no-slip interface condition in (6).Indeed, on the one hand, the presence of the Cauchystress T ´pId of the fluid in ( 4) and ( 5) shows how the velocity and the pressure of the fluid affect the motion of the bodies.On the other hand, considering the Navier-Stokes system (2), ( 3), ( 6) by itself, one can regard the interface condition in (6) as a part of the boundary conditions, which describes the impact of the velocity of each solid on the fluid velocity.The first identity in ( 6) is the standard no-slip boundary condition.

Electromagnetic subsystem
The electromagnetic quantities are determined by the following version of the Maxwell system, together with Ohm's law and the boundary and interface conditions Bptq ¨n " 0 on BΩ, B f ptq ´Bs ptq " 0 on BSptq, Eptq ˆn " 0 on BΩ, ´Ef ptq ´Es ptq ¯ˆn " 0 on BSptq.
In this system we have Ampère's law (7), the Maxwell-Faraday equation ( 8), Gauss's law (9) and Gauss's law for magnetism (10).In comparison to the general situation, these equations have undergone two kinds of reductions, c.f. [27], [28]: First, in the solid region, the system is adjusted to the assumption of the rigid bodies being insulating and second, in the fluid domain, the system is reduced according to the magnetohydrodynamic approximation, see [8,13].A justification for the latter simplification, which is inherent to magnetohydrodynamics, from a physical point of view is for example given in [33,34].The mathematical use of the magnetohydrodynamic approximation consists of the fact that it enables us to summarize the whole electromagnetic problem in Q f into a problem for the magnetic induction B, c.f. the induction equation (29) in our definition of variational solutions below.Once B has been determined, the remaining unknowns E and j are given explicitly through the relations (7) and (11).The quantity J in (7) represents, as in [27], [28], a supplementary external force.
Ohm's law (11) is what links the electromagnetic subsystem to the mechanical subsystem ( 2)-( 6) by describing the influence of the (fluid) velocity u " u f on the electromagnetic quantities in Q f .Since the electrical conductivity σ satisfies σ " σ s " 0 in Q s , it also shows that the electromagnetic fields are not affected by the solid velocities u s i .In the boundary and interface conditions (12), (13) the assumption of the magnetic permeability µ ą 0 being constant in the whole domain Q (c.f.Section 1.2) comes into play.Indeed, while the boundary condition for B in (12) as well as the conditions ( 13) on E are standard, the continuity of B across the interface in (12) is not.However, it is standard to assume continuity of the normal component of B and of the tangential component of 1 µ B and hence, when µ is constant across BSptq, we infer also the second relation in (12).The reason why we cannot allow B to have a jump across the interface is because otherwise it could not be a Sobolev function in Ω, c.f. (24) below.
Both the mechanical subsystem (2)-( 6) and the electromagnetic subsystem ( 7)-( 13) can be studied in their own right, if one considers B and u as a prescribed external forcing term in Ohm's law (11) and the momentum equation (3), respectively.In this article, however, we study the combined system (2)-( 13), in which both B and u are regarded as unknowns, coupled via Ohm's law and the presence of the Lorentz force 1  µ curl B ˆB in the momentum equation.
2 Variational formulation and main result

Notation
The initial positions of the solid bodies are characterized through sets S i 0 Ă Ω, i " 1, ..., m, onto which we impose the conditions Since the motion of the bodies is rigid, we can associate to each body an isometry X i pt, ¨q : R 3 Ñ R 3 , t P r0, T s, such that its position S i ptq at an arbitrary time t P r0, T s can be expressed through the set-valued function In particular, with the notation Xpt; ¨q : S 0 :" :" X i pt; ¨q @i " 1, ..., m, t P r0, T s, the solid region at the time t is given by Sptq :" Xpt; S 0 q.In order to connect the motion of the bodies to the velocity field u we require u to be compatible with the system tS i 0 , X i u m i"1 , i.e. we require the existence of rigid velocity fields u s i pt, ¨q, i " 1, ..., m, such that upt, xq " u s i pt, xq for a.a.t P r0, T s and a.a.x P S i ptq (15) and X i is the unique Carathéodory solution (c.f.[38,Theorem 1.45]) to the initial value problem Finally we denote by T pSq and Y pSq the test function spaces for our variational formulations of the momentum equation and the induction equation, respectively, in Definition 2.1 below.

Weak solutions
We are now in the position to present our variational formulation of the system (2)- (13).With a slight abuse of notation we will write here and in the following sections σ " σ f ą 0, since the quantities containing σ s " 0 are not visible in this weak formulation due to the non-conductivity of the solids.
Definition 2.1 Let T ą 0, let Ω Ă R 3 be a bounded domain and let S 0 " Ť m i"1 S i 0 , where S i 0 Ă Ω for i " 1, ..., m P N satisfy the conditions (14).Assume ν, λ, a, γ, σ, µ P R to satisfy ν, a, σ, µ ą 0, ν `λ ě 0, γ ą Moreover, consider some external forces g, J P L 8 pp0, T q ˆΩq and some initial data 0 ď ρ 0 P L γ pΩq, pρuq 0 P L 1 pΩq, B 0 P L 2 pΩq such that |pρuq 0 | 2 ρ 0 P L 1 pΩq, pρuq 0 " 0 a.e. in tx P Ω : Then the system (2)-( 13) is said to admit a weak solution if there exists a function where each X i pt; ¨q : R 3 Ñ R 3 denotes an isometry, and if there exist functions u P where S " Sp¨q " Xp¨, S 0 q, such that ρ and u, extended by 0 in R 3 zΩ, satisfy the continuity equation and its renormalized form, B t ζpρq `div pζ pρq uq `"ζ 1 pρq ρ ´ζ pρq for any such that the momentum equation and the induction equation, are satisfied for any φ P T pSq and any b P Y pSq, such that the initial conditions ρp0q " ρ 0 , pρuq p0q " pρuq 0 , Bp0q " B 0 , hold true, where the latter two equations are to be understood in the sense that lim for all φ, b P DpΩq with Dpφq " 0 and curl b " 0 in a neighbourhood of S 0 , respectively, and, finally, such that the system tS i 0 , X i u m i"1 is compatible with the velocity field u.
In this definition the compatibility of the velocity field u and the system of rigid bodies leads to some vivid consequences for the solids.First of all, while the bodies are able to touch each other or the domain boundary, the possibility of interpenetrations are ruled out, c.f. [17, Lemma 3.1, Corollary 3.1].Moreover, even though the density does not satisfy a transport equation in the case of a compressible fluid, it still travels along the characteristics of u in the solid part of the domain, c.f. [17,Lemma 3.2].
For definiteness we present these results in the following lemma.
Lemma 2.1 ( [17]) Let Ω Ă R 3 and S i 0 Ă R 3 , i " 1, ..., m P N be bounded domains and let further u P L 2 p0, T ; H 1,2 0 pΩqq be extended by 0 outside of Ω.Moreover, assume u to be compatible with the system tS i 0 , X i u m i"1 where each X i ptq : R 3 Ñ R 3 , t P r0, T s, i " 1, ..., m, denotes an isometry.Then it holds: (i) If, for i ‰ j P t1, ..., mu, there exists τ P r0, T s such that S i pτ q Ş S j pτ q ‰ H then X i ptq " X j ptq for all t P r0, T s.Further, if there exists τ P r0, T s such that S i pτ q Ć Ω, then X i ptq " Id for all t P r0, T s.
(ii) If ρ P L 8 p0, T ; L γ pΩqq, γ ą 1, extended by 0 outside of Ω, satisfies A detailed proof of the assertions (i) and (ii) is given in [17, Lemma 3.1, Corollary 3.1] and [17, Lemma 3.2], respectively.The first part of assertion (i) can be shown directly from the fact that tS i 0 , X i u and tS j 0 , X j u are compatible with the same velocity field u, the second part then follows by regarding R 3 zΩ as a rigid body with the associated rigid velocity field 0. The proof of the assertion (ii) is achieved via a regularization of ρ with respect to the spatial variable and a subsequent application of the regularization method by DiPerna and Lions, c.f. [11], to the continuity equation (32) on compact subsets of the solid time-space domain.

Main result
Our main result yields existence of the weak solutions as introduced in Definition 2.1.
Theorem 2.1 Let T ą 0, assume Ω Ă R 3 to be a simply connected domain of class C 2 Ť C 0,1 and assume S i 0 Ă Ω, i " 1, ..., m P N to be domains of class C 2 Ť C 0,1 which satisfy the conditions (14).Assume moreover the coefficients σ, µ, ν, λ, a, γ P R to satisfy the conditions (19) and the data g, J P L 8 pp0, T q ˆΩq, ρ 0 P L γ pΩq, pρuq 0 P L 1 pΩq and B 0 P L 2 pΩq to satisfy the conditions (20).Then the system (2)-( 13) admits a weak solution tX, ρ, u, Bu in the sense of Definition 2.1 which in addition satisfies the energy inequality for almost all t P r0, T s.
The remainder of the article is devoted to the proof of Theorem 2.1.In the following section we begin with an outline of the main ideas by introducing the approximation method on which the proof is based.

Approximate system
The biggest challenge in the extension of the proof in the incompressible case in [1] to the compressible case lies in the correct construction of the approximate problem.We fix five parameters α, , η ą 0, n P N and ∆t ą 0 and introduce an approximation which consists of five different approximation levels, each of which corresponds to one of the parameters.The approximate system is chosen such that it is easy to solve; a solution to the original system is obtained by passing to the limit in all of the approximation levels.The first three approximation levels, associated to α, and η, respectively, correspond to the approximation used in [17] for the purely mechanical problem: On the α-andlevels, the system is regularized through the addition of an artificial pressure term and multiple further regularization terms.The η-level consists of a penalization method which allows us to pass from a fluid-only system to a system containing both a fluid and rigid bodies.The fourth level, indexed by n, describes a Galerkin method used for solving the approximate momentum equation.Finally, on the fifth level, associated to ∆t, the induction equation is discretized with respect to the time variable while the mechanical part of the problem is split up into a series of time-dependent problems on the small intervals between the discrete times.In the following we present the complete approximate system on the highest approximation level and subsequently give a more explicit description of each included approximation level and its purpose.
Let β ą 0 be sufficiently large such that it satisfies in particular β ą maxt4, γu as in [17,Section 6].Let V n , n P N, denote the Galerkin space spanned by the first n eigenfunctions of the Lamé equation in Ω which constitute an orthonormal basis of L 2 pΩq and an orthogonal basis of H 1,2 0 pΩq, c.f. [37, Lemma 4.33].Then, provided that the approximate system has already been solved up to the (discrete) time pk ´1q∆t for some k " 1, ..., T ∆t , the approximate problem on the interval rpk ´1q∆t, k∆ts consists of finding a solution to the system for all φ P Cprpk ´1q∆t, k∆ts; V n q and b P W k pS ∆t q :" which in addition satisfies the initial conditions ρ ∆t,k ppk ´1q∆t; xq " ρ ∆t,k´1 ppk ´1q∆t; xq, ρ ∆t,1 p0; xq " ρ 0 pxq, x P Ω, ( u ∆t,k ppk ´1q∆t; xq " u ∆t,k´1 ppk ´1q∆t; xq, u ∆t,1 p0; xq " u 0 pxq, x P Ω, (42) Before we proceed with the explanation of the different approximation levels in ( 34)-( 43), let us clarify the notation introduced in this system: For the definition of the set S ∆t pk∆tq in ( 36) and ( 40), we first denote by O i :" pS i 0 q δ :" x P S i 0 : dist `x, BS i 0 ˘ą δ ( the δ-kernel of the initial domain S i 0 of the i-th body, where δ ą 0 is chosen sufficiently small such that for all i " 1, ..., m the δ-neighbourhood pO i q δ " tx P R 3 : distpx, pO i q δ q ă δu of O i coincides with Moreover, we denote by X ∆t,k the unique solution to the initial value problem where R δ ru ∆t,k spt, ¨q :" u ∆t,k pt, ¨q ˚Θδ p¨q and Θ δ denotes a radially symmetric and non-increasing mollifier with respect to the spatial variable.With this notation at hand we define the domain S i ∆t ptq of the i-th approximate solid at an arbitrary time t P rpk ´1q∆t, k∆ts Ă r0, T s by Consequently, the approximate solid region at time t P r0, T s is given by which in particular defines the set S ∆t pk∆tq in (36) and (40).We note that, by construction, S ∆t ptq can be an arbitrary subset of R 3 , while the corresponding approximate solid time-space domain Q s pS ∆t q, defined according to (1), is always a subset of the bounded domain Q.For later use we further remark that S i ∆t ptq, as the δ-neighbourhood of a bounded set, satisfies the cone condition and thus has the property Next, for the definition of the variable viscosity coefficients νpχ k´1 ∆t q and λpχ k´1 ∆t q in the momentum equation (38) we denote the signed distance function of arbitrary sets U Ă R 3 by Further we introduce the signed distance function of the approximate solid area, χ ∆t pt, xq :" db S ∆t ptq pxq, χ k´1 ∆t ptq :" χ ∆t ppk ´1q∆t, ¨q for t P rpk ´1q∆t, k∆ts.

Choosing a convex function H P C 8 pRq such that
Hpzq " 0 for z P p´8, 0s, Hpzq ą 0 for z P p0, `8q, we then define the variable viscosity coefficients by Finally, in the induction equation ( 39) the function ũk´1 ∆t is defined by while the discretized external force J k ∆t is defined by for another mollifier θ ω : R Ñ R and a suitable choice of ω " ωp∆tq, ωp∆tq Ñ 0 for ∆t Ñ 0. We are now in the position to discuss the several approximation levels and the reasons why they are required.We start from the highest level.
The ∆t-level constitutes the level which contains most of the difficulties.It is here where the main novelties of our proof enter, compared to the incompressible setting in [1].The situation presents itself in the following way: On the one hand, the dependence of the test functions ( 17), (18) for the induction equation and the momentum equation on the solution of the system hinders the effort to solve all of the equations in the system simultaneously.While we can deal with the test functions in the momentum equation by means of a penalization method (c.f. the η-level below), the same does not work in case of the induction equation.This suggests to decouple the system by the use of a classical time discretization, e.g.via the Rothe method ([38, Section 8.2]).In this way, at each fixed discrete time we can first determine a velocity field and, from this, the position of the approximate solid.This in turn determines the test functions (40) and solving the discretized induction equation (39) becomes a routine matter.On the other hand, however, the various functions evaluated at different discrete times in a fully discretized system complicate the derivation of a meaningful energy inequality.The author could not find a way to transfer several of the techniques known for the continuous compressible Navier-Stokes system (c.f.[37, Sections 7.6.5,7.6.6,7.7.4.2]) -in particular, the proof of the non-negativity of the density -to the discrete case and it did not seem to be possible to derive the uniform bounds required for the limit passage with respect to ∆t Ñ 0. Our solution to this dilemma consists of considering, instead of a strictly discretized system, a hybrid system in which the induction equation ( 39) is indeed discretized by the Rothe method, while the continuity equation and the momentum equation are solved as continuous equations on the small intervals between each pair of consecutive discrete times, c.f (37) and (38).Through this, the solution dependence of the test functions in the induction equation can be handled as in the fully discrete system, while the mechanical part of the energy inequality -with the density bounded away from zero -can be derived as in the strictly continuous case.Moreover, under the consideration of piecewise linear interpolants of the discrete functions, the discrete induction equation ( 39) also leads to a continuous energy estimate, which can be combined with the mechanical estimate to obtain the full energy inequality, c.f. Section 5.1.A hybrid approximation scheme was already used in our proof in the incompressible case [1].In that case, however, the major part of the system could be discretized in time while only the transport equation for the characteristic function of the solid domain had to be treated as a continuous problem on small time intervals.The idea for the latter procedure, in turn, stems from [26].
The Galerkin method carried out on the n-level is used to solve the continuous momentum equation (37) on the small time intervals from the ∆t-level by a standard procedure.The Galerkin-regularity of the velocity field furthermore helps us during the limit passage with respect to ∆t Ñ 0, c.f. (85) below.
After letting n tend to 8 we find ourselves in the same situation as in the approximation of the exclusively mechanical system in [17, Section 6].Indeed, the remaining three approximation levels correspond directly to the three level approximation scheme used in that article.Hence, for the mechanical part of our problem we can follow exactly the strategy used therein.Moreover, the limit passages in the induction equation from here on do not contain any new difficulties anymore.Consequently, after the limit passage in n the rest of the proof will become a routine matter.The penalization method on the η-level -c.f.Section 7 -is the same as the one used for the fluidrigid bodies system in [17] and was, before that, also used for example for the corresponding two dimensional problem in [40].The idea behind it is to approximate the entirety of the fluid and the rigid bodies by a fluid in the whole domain with viscosity tending to infinity in the later solid regions.Mathematically this is implemented through the variable viscosity coefficients (49).Due to the choice of the function H in (48) these coefficients blow up in the approximate solid region once we let η tend to 0 and, thanks to the energy inequality, this will cause the limit velocity field u to coincide with a rigid velocity field in each body.Moreover, the positions S i ptq of the bodies in the η-limit are determined through the flow curves of R δ rus, c.f. ( 44) and ( 46).This regularized velocity field has the useful property that, for any domain c.f. [17, Remark 6.1].Hence R δ rus coincides with u itself in the sets O i ptq Ă S i ptq, in which Dpuq " 0. Consequently, the rigid velocity fields coinciding with u in S i ptq also coincide, in O i ptq, with the velocity field R δ rus which determines the motion of the bodies.In particular, this shows that the bodies S i ptq are indeed rigid.On the -level, the continuity equation ( 37) is regularized through the additional Laplacian ∆ρ ∆t,k .
The additional quantity p∇u ∆t,k ∇ρ ∆t,k q in the momentum equation (38)  ˇˇ2 curl B k ∆t q in the induction equation (39) fulfills, as in the incompressible setting in [1], the same purpose for the magnetic part of the mixed terms so that we are able to derive uniform bounds from the energy inequality nevertheless, c.f. Section 5.1.We remark that this control of the mixed terms was also the motivation for the definition of the velocity field (50) in the discrete induction equation: Indeed, defining this quantity as a mean value of the velocity field obtained from the momentum equation on the intervals rpk´2q∆t, pk´1q∆ts, we can absorb it into the left-hand side of the energy inequality thanks to the above-mentioned regularization terms.If instead the term was defined, more intuitively, as a pointwise evaluation of u ∆t,k , we would not be able to handle it.The last regularization term in (39), the 4-th curl of B k ∆t , enables us to express the induction equation via some weakly continuous operator on Y k pS ∆t q.Seeing that this operator is moreover coercive, we will then be able to infer the existence of B k ∆t , c.f. Section 4.2.Finally, on the α-level, the artificial pressure term αρ β ∆t,k is added to the momentum equation (38).Again this method is already well-known from the general existence theory for the compressible Navier-Stokes system, c.f. [37, Section 7.3.8].The artificial pressure gives us an additional amount of integrability of the density and its gradient, required to pass to the limit in the term p∇u ∆t,k ∇ρ ∆t,k q from the -level, c.f. [37, Section 7.8.2].It furthermore simplifies the limit passage with respect to Ñ 0, since the additional integrability allows for the use of the regularization technique by DiPerna and Lions, c.f. [37, Lemma 6.8, Lemma 6.9].

Existence of the approximate solution
We begin the proof of Theorem 2.1 by showing the existence of a solution to the approximate problem ( 34)-( 43) on the highest approximation level.

Existence of the density and velocity
The existence of the density and the velocity field on the Galerkin level can be shown by classical methods, c.f. for example [37,Section 7.7].More precisely, the continuity equation (37) and the momentum equation ( 38) can be solved simultaneously by means of a fixed point argument: For fixed w P Cprpk ´1q∆t, k∆ts; V n q we consider the Neumann problem ∇ρ ¨n| BΩ " 0, ρppk ´1q∆t, ¨q " ρ ∆t,k´1 ppk ´1q∆t; ¨q in Ω, (54) It is well known (c.f.Lemma 10.1 in the appendix) that ( 53)-( 55) admits a unique solution ρ " ρpwq P C ´rpk ´1q∆t, k∆ts; C ´ p∇u∇ρpwqq ¨φ ´ |w| 2 u ¨φ dx in rpk ´1q∆t, k∆ts, uppk ´1q∆t, ¨q "u ∆t,k´1 ppk ´1q∆t; ¨q in Ω, for all φ P Cprpk ´1q∆t, k∆ts; V n q.Under exploitation of the fact that, by (56), ρpwq is bounded away from 0 and the linearity of the problem, it follows from classical methods that this problem admits a unique solution u " upwq P Cprpk ´1q∆t, k∆ts; V n q.We can thus define an operator T : C prpk ´1q∆t, k∆ts; V n q Ñ C prpk ´1q∆t, k∆ts; V n q , Tpwq :" u.
It is easy to see that T is continuous and compact and, by an energy estimate, fixed points of sT for s P r0, 1s are bounded in Cprpk ´1q∆t, k∆ts; V n q, uniformly with respect to s.Under these conditions the Schaefer fixed point theorem (see [14, Section 9.2.2,Theorem 4]) tells us that T possesses a fixed point u ∆t,k P Cprpk ´1q∆t, k∆ts; V n q, which constitutes the desired solution to the initial value problem (38), (42).Furthermore, by construction, the associated density ρ ∆t,k :" ρpu ∆t,k q is the desired solution to the corresponding initial value problem (37), (41) for the density.

Existence of the magnetic induction
The existence of the magnetic induction is obtained as in the incompressible case, c.f. [1, Section 3].We equip the space Y k pS ∆t q with the norm } ¨}H 2,2 pΩq and express the identity (39) through the equation where the operator A : Y k pS ∆t q Ñ pY k pS ∆t qq ˚and the right-hand side f P pY k pS ∆t qq ˚are defined by xApBq, by pY k pS ∆t qq ˚ˆY k pS ∆t q :" for any B, b P Y k pS ∆t q.The operator A is weakly continuous and coercive and consequently surjective from Y k pS ∆t q onto pY k pS ∆t qq ˚, c.f. [22, Theorem 1.2].In particular, there exists a function B k ∆t P Y k pS ∆t q which satisfies (57) and thus the induction equation (39) for all b P Y k pS ∆t q.Finally, by the Helmholtz decomposition, see [43,Theorem 4.2], we infer that (39) does not only hold for b P Y k pS ∆t q but also for the (not divergence-free) test functions b P W k pS ∆t q.Altogether, we have shown the following result: Proposition 4.1 Let n P N, ∆t, η, , α ą 0 such that T ∆t P N and let β ą maxt4, γu be sufficiently large.Assume the conditions of Theorem 2.1 to be satisfied.Moreover, assume that where P n denotes the orthogonal projection of L 2 pΩq onto V n .Finally, for all k " 1, ..., T ∆t , let J k ∆t be defined by (51).Then, for each k " 1, ..., T ∆t , there exist functions which satisfy the continuity equation (37), the momentum equation (38) for all test functions φ P C prpk ´1q∆t, k∆ts; V n q and the induction equation (39) for all test functions b P W k pS ∆t q as well as the initial conditions (41)- (43).

Limit passage in the time discretization
We continue by passing to the limit with respect to ∆t Ñ 0. As in [1], we first need to assemble the functions constructed in Section 4, defined up to now only on small time intervals or in discrete time points, to functions defined on the whole time interval r0, T s.More precisely, for functions f ∆t,k , defined on rpk ´1q∆t, k∆ts ˆΩ for k " 1, ..., T ∆t , we introduce the assembled functions  37), the momentum equation (38), and the initial conditions ( 41)- (43), it follows from the definition of ρ ∆t and u ∆t in (58) as well as of B ∆t , B ∆t and B 1 ∆t in (59)-(61) that these functions solve the continuity equation B t ρ ∆t `div pρ ∆t u ∆t q " ρ ∆t a.e. in r0, T s ˆΩ, the momentum equation for any φ P Cpr0, T s; V n q and the initial conditions Furthermore, from X ∆t,k being the unique solution to the initial value problem ( 44), (45), it follows that X ∆t is the unique solution to Finally, we consider functions b P L 4 ´0, T ; H 2,2 0 pΩq ¯such that bpτ q P W l pS ∆t q for a.a.τ P rpl ´1q∆t, l∆ts, l " 1, ..., and realize that, after a density argument, the discrete induction equation (39) at time k∆t can be tested by bptq for almost all t P rpk´1q∆t, k∆ts.After integration over rpk´1q∆t, k∆ts and summation over k this yields ´ curl `curl B ∆t ˘¨curl pcurl bq dxdt. (67)

Energy inequality on the ∆t-level
In contrast to the incompressible setting in [1] we have to combine the discrete induction equation (39) with the continuous Navier-Stokes equations ( 37), (38) in a suitable way in order to derive an energy inequality at the ∆t-level.We pick an arbitrary t P p0, T s and choose k P 1, ..., T ∆t ( , ξ P r0, ∆tq such that t " k∆t ´ξ.For the magnetic part of the energy inequality we test the induction equation ( 39) by 1 µ B k ∆t , which leads to Since corresponding estimates hold true also for all time indices l " 1, ..., k ´1 we can integrate (discretely) over the interval r0, ts, which yields under exploitation of Hölder's and Young's inequalities in the last estimate.On the right-hand side of (69) we further estimate, due to the definition of ũl´1 ∆t in (50) and Jensen's inequality, |u ∆t,l´1 pτ q| 4 dτ dx.
Moreover, a direct calculation yields Hence, absorbing the } curl B l ∆t } 4 L 4 pΩq -terms in (69) into the left-hand side and expressing the sums as integrals, we end up with For the mechanical part of the energy inequality we test the continuity equation (62) by 1 2 |u ∆t | 2 , the momentum equation ( 63) by u ∆t , add up the resulting equations and obtain where the last estimate uses Hölder's inequality, Young's inequality and the same inequality (70) as in the estimate of the corresponding term in the induction equation.Adding the inequality (72) to the inequality (71) and absorbing multiple terms from the right-hand side into the left-hand side, we finally get the energy inequality where the constant c ą 0 is independent of ∆t and t.In particular, by use of the Gronwall Lemma and the estimates for the solution to the Neumann problem for the density, c.f. Lemma 10.1, we find a constant c ą 0, independent of ∆t, such that the following bounds hold true: The bounds for the magnetic induction in L 8 p0, T ; L 2 pΩqq in (75)-(77) follow from the choice t " k∆t, k " 1, ..., T ∆t in the energy inequality (73), for which it holds B ∆t,}¨} ptq " }B k ∆t } 2 L 2 pΩq .For a bound of the time derivative of u ∆t we introduce the operator and denote N pρ ∆t , u ∆t q :" ´div pρ ∆t u ∆t b u ∆t q ´a∇ρ γ ∆t ´α∇ρ β ∆t `div ´2ν ´χk´1 ∆t ¯D pu ∆t q ∇ ´λ ´χk´1 ∆t ¯div pu ∆t q ¯`ρ ∆t g ` This together with the uniform bound of ρ ∆t away from 0 and the uniform bounds (74)-(77) leads to the estimate For more details on the derivation of (78) we refer to [37, Section 7.7] and in particular [37, Section 7.7.4.1].The bounds (74)-( 77) and ( 78) and the Aubin-Lions Lemma imply the existence of functions such that, after the extraction of a subsequence, it holds Here, the fact that the weak limits of the different interpolants coincide follows from Lemma 10.2.Moreover the boundary conditions of the limit functions in (79) and (80) follow directly from the corresponding boundary conditions on the ∆t-level, c.f. Proposition 4.1 and the definition of the space Y k pS ∆t q in (36).Furthermore, the external force J ∆t , discretized via (51), converges to its original time-dependent counterpart, J ∆t Ñ J in L q pQq @1 ď q ă 8.
Finally, X ∆t , as the solution to the initial value problem (65), satisfies the conditions of Lemma 10.3 which tells us that where Sptq :" pXpt; Oqq δ and X represents the solution to dXpt; xq dt " R δ rus pt, Xpt; xqq , for t P r0, T s, Xp0; xq " x, for x P R 3 .

Continuity equation
Due to the convergences (81)-(83) of the density and the velocity we can pass to the limit in the continuity equation (62) and infer that the limit functions ρ and u solve the initial value problem B t ρ `∇ pρuq " ∆ρ a.e. in Q, ρp0q " ρ 0 .

Induction equation
We first show convergence of the quantity ũ1 ∆t from the mixed term in the discrete induction equation (67).We fix an arbitrary point pt, xq P Q and, for each sufficiently small ∆t ą 0, we choose k ∆t P t2, ..., T   due to the uniform convergence (83) of u ∆t .Moreover, the uniform bound of u ∆t in (74) shows equiintegrability of |ũ 1 ∆t ´u∆t | q for any 1 ď q ă 8.This together with the pointwise convergence (85) gives us the conditions for the Vitali convergence theorem and we infer that Further, due to the uniform bounds (76), we find a function z P L 4 3 pQq such that, possibly after the extraction of a suitable subsequence, it holds We remark that, since the quantity z will vanish from the system as soon as we let tend to zero, there is no need to specify the form of the limit function z in (87).Now we test (67) by an arbitrary function b P Y pSq.This is possible, since b is curl-free in an open neighbourhood of Q s pSq and so, by the uniform convergence (84) of the signed distance function, it also satisfies curl bptq " 0 a.e. in S ∆t pk∆tq for a.a.t P rpk ´1q∆t, k∆ts, k " 1, ..., T ∆t (88) for all sufficiently small ∆t ą 0. Hence, for all such t, k and ∆t it holds bptq P W k pS ∆t q and thus b is indeed an admissible test function in (67), c.f. (66).Using the strong convergence (86) of ũ1 ∆t , we can pass to the limit with respect to ∆t Ñ 0 in the resulting equation and obtain the induction equation in the limit, ´ curl pcurl Bq ¨curl pcurl bq dxdt @b P Y pSq.

Momentum equation
In order to pass to the limit in the momentum equation, it remains to show convergence of the piecewise constant Lorentz force.This is achieved by the same arguments as in the incompressible case, c.f. [1, Section 4.2]: The uniform bounds ( 76) and (77) allow us to extract suitable subsequences and find functions z 1 , z 2 P L 2 pQq with the properties that curl B 1 ∆t Our goal is to identify the limit functions z 1 and z 2 as Since, according to (47), it holds it is sufficient to do so in Q s pSq and Q f pSq.Focussing at first on the fluid part of the domain, we consider any c, d P r0, T s and any ball U Ă Ω such that I ˆU :" pc, dq ˆU Ă Q f pSq.In particular, any function b P L 4 pI; H 2,2 0 pU qq, extended by 0 outside of I ˆU , is curl-free in an open neighbourhood of Q s pSq.Hence, the uniform convergence (84) of the signed distance function implies that, for sufficiently small ∆t ą 0, any such b also satisfies the curl-free condition (88) in S ∆t pk∆tq, k " 1, ..., T ∆t , and therefore constitutes an admissible test function in the induction equation (67).Using it as such we deduce the uniform bound Since, by (77), B 1 ∆t is bounded uniformly in L 2 pQq, we can thus apply the discrete Aubin-Lions Lemma [12, Theorem 1] and conclude that This strong convergence implies that z 1 " z 2 " curl B ˆB a.e. in I ˆU and hence in Q f pSq. (91) In order to show the equation ( 90) also in the solid region, we consider another arbitrary pair of an interval I Ă p0, T q and a ball U Ă Ω, this time satisfying I ˆU Ă Q s pSq.From the fact that curl B k ∆t " 0 in S k ∆t Ş Ω for k " 1, ..., T ∆t , and the uniform convergence (84) of the signed distance function it follows that, for any sufficiently small ∆t ą 0, curl B 1 ∆t " curl B ∆t " 0 a.e. in I ˆU.
Thus, letting ∆t tend to zero, we infer that z 1 " z 2 " 0 " curl B " curl B ˆB a.e. in I ˆU and hence in Q s pSq. (92) In combination with the corresponding equality (91) in Q f pSq, we have therefore shown the desired identification (90) of z 1 and z 2 .We remark that (92) moreover shows that the solid region Q s pSq in the limit is again insulating.Next we exploit the uniform convergence (84) of the signed distance function together with the definition (49) of the variable viscosity coefficients as smooth functions of χ 1 ∆t to infer that This, in combination with the convergence of the Lorentz force, c.f. ( 89), (90) and the convergences (81)-(83) of ρ ∆t and u ∆t , allows us to pass to the limit in the momentum equation ( 63) and infer the equation Ω pρu b uq : Dpφq `´aρ γ `αρ β ¯div φ ´2ν pχq D puq : Dpφq ´λ pχq div puq div φ `ρg ¨φ `1 µ pcurlB ˆBq ¨φ ´ |u| 2 u ¨φ ´ p∇u∇ρq ¨φ dxdt @φ P C pr0, T s; V n q .

Limit passage in the energy inequality
We choose an arbitrary t P p0, T s and k P t1, ..., T ∆t u, such that t " k∆t ´ξ for some ξ P r0, ∆tq.Subsequently, we sum the inequality (68) over l " 1, ..., k and add the first inequality in (72) to obtain Here, on the left-hand side we drop the term ρ γ´2 ∆t |∇ρ ∆t | 2 and pass to the limit by exploiting, in particular, the weak lower semicontinuity of norms and the strong convergence (93) of the variable viscosity coefficients.Moreover, on the right-hand side of (94) we can carry out the limit passage by exploiting the convergence (89), (90) of the Lorentz force.Altogether we obtain Hence we have proved Proposition 5.1 Let n P N, η, , α ą 0, β ą maxt4, γu sufficiently large and let the assumptions of Theorem 2.1 be satisfied.Furthermore, let Under these conditions there exist functions X n : R 3 Ñ R 3 , z n P L 4 3 pp0, T q ˆΩq and for S n " S n p¨q " pX n p¨; Oqq δ , which satisfy dX n pt; xq dt " R δ ru n s pt, X n pt; xqq , for t P r0, T s, X n p0; xq " x, for x P R 3 , B t ρ n `div pρ n u n q " ∆ρ n a.e. in Q (97) ´ p∇u n ∇ρ n q ¨φ dxdt, (98) ´ curl pcurl B n q : curl pcurl bq dxdt, where χ n pt, xq :" db Snptq pxq, for any φ P Cpr0, T s; V n q and any b P Y pS n q.Further, these functions satisfy the initial conditions of which the latter identity can be understood in the sense of (31), as well as the energy inequality for almost all t P r0, T s.

Limit passage in the Galerkin method
Next, we pass to the limit with respect to n Ñ 8. Using Lebesgue interpolation, we infer from the energy inequality (100) the existence of a constant cp , αq ą 0, independent of n, such that L 8 p0,T ;L β pΩqq }ρ n } 3 5 L β p0,T ;L 3β pΩqq ď cp , αq.
Moreover, from the classical L p -L q regularity results for parabolic equations, c.f. [37, Lemma 7.37, Lemma 7.38, Section 7.8.2],we infer that ρ n as the solution to the regularized continuity equation (97) satisfies the estimates and a constant c ą 0 independent of n, η and .The uniform bounds (101),(102) and the energy inequality (100), together with the Aubin-Lions Lemma, allow us to extract suitable subsequences and find functions z P L 4 3 pQq and 0 ď ρ P " ψ P L 8 ´0, T ; L β pΩq ¯č L r `0, T ; W 1,r pΩq ˘č L r `0, T ; W 2,r pΩq ˘: with the properties that The boundary conditions of the limit functions in (104)-( 106) follow directly from the corresponding boundary conditions (95) and (96) of ∇ρ n and B n and the fact that u n P V n vanishes on BΩ for all n P N. Finally, the initial value problem (65), solved by X n , yields that the conditions of Lemma 10.3 are satisfied.Hence where Sptq :" pXpt; Oqq δ and X denotes the unique solution to the initial value problem dXpt; xq dt " R δ rus pt, Xpt; xqq , for t P r0, T s, Xp0; xq " x, for x P R 3 .

Continuity equation
The convergences (107)-(109) of ρ n and u n allow us to pass to the limit in the continuity equation (97).Consequently, the limit functions ρ and u satisfy the continuity equation B t ρ `div pρuq " ∆ρ a.e. in Q, ρp0q " ρ 0 .

Induction equation
At this stage -as well as in the later sections -the limit passage in the induction equation does not differ from the incompressible case, c.f. [1].For the convenience of the reader, we present the arguments here: We begin by making sure that test functions b P Y pS, Xq for the limit equation are also admissible on the n-level.To this end we fix some arbitrary ω ą 0. Then the uniform convergence (110) of the signed distance function implies the existence of N " N pωq ą 0, such that for any function b P Y pSq with curl b " 0 in an ω-neighbourhood of Q s pSq it also holds curl b " 0 in an ω 2 -neighbourhood of Q s pS n q and thus b P Y pS n q @n ě N.
In particular, for any interval I Ă p0, T q and any ball U Ă Ω, such that I ˆU Ă Q f pSq there exists N " N pI ˆU q ą 0 such that where b has been extended by 0 outside of I ˆU .Next, an interpolation between L 8 p0, T ; L 2 pΩqq and L 2 p0, T ; L 6 pΩqq provides a uniform bound of B n in L 3 pQq.In combination with the bounds of u n and curl B n in L 2 pQq this yields the existence of functions z 3 , z 4 P L With the aim of identifying z 3 and z 4 we pick an arbitrary interval I Ă p0, T q and an arbitrary ball U Ă Ω with the property I ˆU Ă Q f pSq.From (111) we know that for any sufficiently large n P N the induction equation (99) may be tested by all functions of the form ψb, where ψ P DpIq, b P DpU q.Doing so, we obtain the dual estimate with a constant c ą 0 depending on b but not on n.This, together with the Arzelá-Ascoli theorem, yields Combining this with the weak convergences of u n and B n in L 2 p0, T ; H 1,2 pΩqq we realize that Due to the uniform convergence (110) of the signed distance function and the fact that curl B n " 0 in Q s pS n q, the former one of these identities also holds true in Q s pSq, Since moreover test functions b P Y pSq are curl-free in Q s pSq, we end up with for any b P Y pSq.The relations (112) and ( 114) allow us to pass to the limit in the mixed term of the induction equations and hence, letting n tend to infinity in (99), we see that ´ curl pcurl Bq ¨curl pcurl bq dxdt for some small t 0 " t 0 pbq ą 0. This yields the initial condition Bp0q " B 0 in the sense of (31).

Momentum equation
By the same methods as for the (purely mechanical) compressible Navier-Stokes system, c.f. [37, Section 7.8.2],we derive strong convergence of the momentum function in the Galerkin limit: Recalling that P n denotes the orthogonal projection of L 2 pΩq onto V n , we test the momentum equation ( -after a density argument -by P n pφq for an arbitrary function φ P L 5β´3 β´3 p0, T ; H 2,2 0 pΩqq.Since in particular , this results in the dual estimate Consequently, from the Aubin-Lions Lemma, we conclude that Therefore, as u n converges weakly in L 2 p0, T ; H 1,2 pΩqq, we infer that, for example, Moreover, the bound of u n in L 4 pQq implies the existence of some z P L 4 3 pQq such that, for a chosen subsequence, it holds Using further the uniform convergence (110) of the signed distance function for passing to the limit in the variable viscosity coefficients and the relations (112), (114) which allow us to pass to the limit in the Lorentz force, we can now let n tend to infinity in the momentum equation ( 98 for any φ P C 1 0 pr0, T s; V N q with fixed N P N. Since 0 pΩq, we finally conclude that (115) also holds true for any φ P DpQq.Using further the weak lower semicontinuity of norms to pass to the limit in both the energy inequality (100) and the uniform bounds (102) we have proved the following proposition: Proposition 6.1 Let η, , α ą 0, β ą maxt4, γu sufficiently large and let the assumptions of Theorem 2.1 be satisfied.Furthermore, let Under these conditions there exist functions X η : R 3 Ñ R 3 , u η P L 2 p0, T ; H 1,2 0 pΩqq, z η , zη P L 4 3 pQq and 0 ď ρ η P " ψ P L 8 ´0, T ; L β pΩq ¯č L r `0, T ; W 1,r pΩq ˘č L r `0, T ; W 2,r pΩq ˘: for r ą 2, r ą 1 as in (103) and S η " S η p¨q " pX η p¨; Oqq δ , which satisfy dX η pt; xq dt " R δ ru η s pt, X η pt; xqq , for t P r0, T s, X η p0; xq " x, for x P R 3 , and Dpφq `´aρ γ η `αρ β η ¯div φ ´2ν pχ η q Dpu η q : Dpφq ´λ pχ η q div u η div φ `ρη g ¨φ `1 µ pcurl B η ˆBη q ¨φ ´ zη ¨φ ´ p∇u η ∇ρ η q ¨φ dxdt, ´ curl pcurl B η q : curl pcurl bq dxdt, where χ η pt, xq :" db Sηptq pxq, for any φ P DpQq and any b P Y pS η q.Further, these functions satisfy the initial conditions ρ η p0q " ρ 0 , pρ η u η q p0q " pρuq 0 , B η p0q " B 0 , of which the latter identity can be understood in the sense of (31), as well as the energy inequality for almost all t P r0, T s and the estimate with a constant c ą 0 independent of η and .
From this point on, the remainder of the proof of the main result is straight forward: In the mechanical part of the problem we can follow precisely the arguments from [17,, the additional Lorentz force (c.f.[41]) and regularization term in the momentum equation do not cause any essential further difficulties.In the induction equation, each limit passage from now on can be carried out as in the incompressible case in [1] and thus essentially as in Section 6.2.However, for the convenience of the reader, we will sketch the main arguments for the remaining three limit passages in the following sections.

Limit passage in the penalization method
We continue by passing to the limit with respect to η Ñ 0. Exactly as in the limit passage with respect to n Ñ 8 in Section 6 we can, due to the energy inequality (121) and the uniform bound (122), extract suitable subsequences and find functions z, z P L 0 ď ρ P " ψ P L 8 ´0, T ; L β pΩq ¯č L r `0, T ; W 1,r pΩq ˘č L r `0, T ; W 2,r pΩq ˘: such that c.f. Lemma 10.3.In particular, this lemma also implies that χ η Ñ χ :" db Sp¨q p¨q in C `r0, T s; C loc `R3 ˘˘.

Continuity equation
Making use of the convergences (126)-(128) of ρ η and u η , we can pass to the limit in the continuity equation ( 119) and ensure that B t ρ `div pρuq " ∆ρ a.e. in Q, ρp0q " ρ 0 .
Moreover, this pointwise identity can be renormalized by multiplying it by ζ 1 pρ q for an arbitrary convex function ζ P C 2 pr0, `8qq.Since ζ 2 ě 0, this yields almost everywhere in Q.This relation will turn out useful in the limit passage with respect to Ñ 0 in Section 8.

Induction equation
For the limit passage in the induction equation we can argue exactly as in the limit passage with respect to n Ñ 8 in Section 6.2 to show strong convergence of B η in the fluid domain.Hence, we can pass to the limit in (120) and obtain the identity ´ curl pcurl Bq ¨curl pcurl bq dxdt @b P Y pSq.
Moreover, the initial condition Bp0q " B 0 also follows as in Section 6.2.

Momentum equation and compatibility of the velocity field
From the uniform bounds given by the energy inequality (121) we further infer the existence of z 5 P L 6 5 pQq such that, for a chosen subsequence, ρ η u η b η á z 5 in L 6 5 pQq.
For the limit passage in the momentum equation we need to identify z 5 : The choice of test functions φ P T pSq, which satisfy Dpφq " 0 in a neighbourhood of Q s pSq, allows us to control the variable viscosity coefficients νpχ η q and λpχ η q in the momentum equation (120), since these remain bounded in the fluid region according to their definition in (48), (49).This enables us to deduce strong convergence of the momentum function ρ η u η in the fluid domain similarly to the strong convergence (113) of the magnetic induction in the Galerkin limit.Indeed, we fix an arbitrary interval I Ă p0, T q and an arbitrary ball U Ă Ω such that I ˆU Ă Q f pSq and deduce from the momentum equation, for any Φ P DpU q, the dual estimate for a constant c ą 0 depending on Φ but not on η.From the Arzelá-Ascoli theorem it follows that ρ η u η Ñ ρu in C weak ´I; L 2β β`1 pU q ¯and hence in L 2 `I; H ´1,2 pU q ˘, which implies that Moreover, since νpχ η q and λpχ η q blow up in the solid part of the domain, the energy inequality (121) shows that Dpuq " 0 a.e. in Q s pSq.
Hence, there are rigid velocity fields u s i which coincide with u almost everywhere in the δ-neighbourhoods S i ptq :" pXpt; O i qq δ of the sets Xpt; O i q.Consequently, due to the property (52) of the regularized velocity field R δ rus, we can replace R δ rus in the initial value problem (129) by u s i for x P O i .The combination of the latter two conditions at first yields compatibility (c.f. ( 15), ( 16)) of u with the system tO i , X i u m i"1 , where each X i ptq denotes an isometry which coincides with Xptq in O i .However, the fact that each X i ptq is an isometry implies that S i ptq " `Xi `t; O i ˘˘δ " X i `t; S i 0 ˘and thus upt, ¨q " u s i pt, ¨q a.e. in X i `t; S i ˘, Consequently, we infer that u is even compatible with tS i 0 , X i u m i"1 .Finally, the initial condition pρuqp0q " pρuq 0 , in the sense of (31), follows by similar arguments as the initial condition for the magnetic induction in Section 6.2.

Energy inequality
We drop, among other non-negative terms, the variable parts of the viscosity coefficients on the lefthand side of the energy inequality (121) and pass to the limit to see that for almost all t P r0, T s.

Limit passage in the regularization terms
The next step is the limit passage with respect to Ñ 0. The energy inequality (134) yields the bounds needed for Corollary 10.1, which implies the existence of isometries X i ptq : R 3 Ñ R 3 such that We write X : r0, T s ˆS0 Ñ R 3 , Xptq| S i 0 :" X i ptq and S " Sp¨q :" Xp¨; S 0 q.Further, testing the continuity equation (130) by ρ , we see that 1 2 }∇ρ } L 2 pQq ď c for a constant c ą 0 independent of .This, together with the energy inequality (134), yields the existence of functions 0 ď ρ P L 8 p0, T ; L β pΩqq and such that for certain extracted subsequences it holds The boundary conditions of the limit functions in (135) and (136) follow directly from the corresponding boundary conditions of the velocity field and the magnetic induction in ( 123) and ( 125) on the -level.Moreover, the velocity field u is compatible with the system tS i 0 , X i u m i"1 , c.f. Corollary 10.1.

Continuity equation
Similarly to the strong convergence (113) of the magnetic induction, we deduce from the continuity equation ( 130) that ρ even converges to ρ in C weak pr0, T s; L β pΩqq.This, together with the vanishing artificial viscosity term, c.f. (137), is sufficient to pass to the limit in (130) and to obtain In fact, as at this stage of the approximation it holds ρ P L 2 pQq, u P L 2 p0, T ; H 1,2 0 pΩqq, we can use the regularization procedure by DiPerna and Lions, c.f. [37, Lemma 6.8, Lemma 6.9], to see that ρ and u, extended by 0 outside of Ω, even satisfy the renormalized continuity equation ( 26), (27).This in turn implies that ρ P Cpr0, T s; L 1 pΩqq, c.f. [37, Lemma 6.15], and ρ satisfies the initial condition ρp0q " ρ 0 .

Induction equation
The regularization terms in the induction equation vanish as tends to 0 according to the convergences (137).Apart from that we can argue exactly as in the Galerkin limit in Section 6.2 to pass to the limit with respect to Ñ 0 in (132) and to infer that

Momentum equation
In order to pass to the limit in the pressure terms, we first consider an arbitrary compact set K Ă Q f pSq.Denoting by B Ω the Bogovskii operator in Ω (c.f.[37, Section 3.
where Φ P DpQ f pSqq is a cut-off function equal to 1 in K.This procedure leads to a bound of ρ in L β`1 pKq uniformly in , c.f. [17, Lemma 8.1] and the references therein.These bounds in turn allow us to find With the aim of identifying these limit functions we set, for arbitrary Φ P DpQ f pSqq, where ∆ ´1 denotes the inverse Laplacian on R 3 , c.f. [19,Section 10.16].We compare the momentum equation ( 133) on the -level, tested by φ , to a corresponding limit identity, tested by φ.This enables us to deduce the effective viscous flux identity for all Φ P DpQ f pSqq, c.f. [17, Lemma 8.2] and the references therein.Moreover, after a density argument, we can consider the choice ζpsq " s lnpsq in both the renormalized continuity equations (131) on the -level and (26) in the limit.A comparison between the resulting identities then leads us to for τ P r0, T s, where the last inequality follows from the effective viscous flux identity and the monotonicity of the mapping s Þ Ñ as γ `αs β , as well as the fact that div u " 0 in Q s pSq.Due to the strict convexity of ζ this estimate implies pointwise convergence of ρ in p0, T q ˆΩ, c.f. [19, Theorem 10.20], and hence z 6 " ρ γ , z 7 " ρ β almost everywhere in Q f pSq.In the remaining terms of the momentum equation ( 133) we can pass to the limit as during the past limit passages.We end up with

Energy inequality
Neglecting the regularization terms on the left-hand side of the energy inequality (134) on the -level, we can pass to the limit with respect to Ñ 0 and obtain for almost every t P r0, T s.

Limit passage in the artificial pressure
Finally it remains to pass to the limit with respect to α Ñ 0. We now consider initial data ρ 0 , pρuq 0 and B 0 as in Theorem 2.1 and construct -c.f.[21, Section 4] -the initial data ρ 0,α , pρuq 0,α and B 0,α on the α-level (c.f. ( 116), ( 117)) in such a way that Since the energy inequality (134) provides the conditions for Corollary 10.1, we obtain the existence of isometries X i ptq : R 3 Ñ R 3 such that We set and S " Sp¨q :" Xp¨; S 0 q.Then from the energy inequality (144) we obtain the existence of functions 0 ď ρ PL 8 p0, T ; L γ pΩqq, such that for suitable subsequences it holds The boundary conditions of the limit functions in (147) and (148) follow directly from the corresponding boundary conditions of the velocity field and the magnetic induction in ( 135) and ( 136) on the α-level.Moreover, again due to Corollary 10.1, is compatible with the family tS i 0 , X i u m i"1 . (149)

Continuity equation
After using the continuity equation ( 138) to deduce convergence of ρ α in C weak pr0, T s; L γ pΩqq, we pass to the limit in (138) and obtain The proof of the renormalized continuity equation however needs to be postponed to Section 9.3 below, since at this stage ρ does not have the L 2 pQq-regularity required for the regularization technique by DiPerna and Lions anymore.

Induction equation
For the limit passage in the induction equation (139) we argue exactly as in Section 6.2 and end up with

Momentum equation
For the limit passage in the pressure terms the strategy used during the limit passage with respect to Ñ 0 in Section 8.3 needs to be modified to make up for the lower integrability of the density; the main ideas however remain the same.First we test the momentum equation (143) by functions of the form (140) with the density replaced by (a cut-off and smoothened version of) ρ θ α , θ ą 0. Choosing θ ą 0 sufficiently small, we find that, for any compact K Ă Q f pSq, ρ α and α where T k pρq denotes a weak L 1 -limit of T k pρ α q, holds true for any φ P DpQ f pSqq.
The effective viscous flux identity (152) can be shown by a comparison between the momentum equation (143) on the α-level and a corresponding limit identity, tested by suitably modified variants of the functions φ and φ in (141) with the density replaced by T k pρ α q and T k pρq, respectively.The details, in the case without rigid bodies, are given e.g. in [21,Section 4.3], the adjustment to the fluid-structure case poses no further difficulties.The proof of the bound (153) of the oscillation defect measure is split up into an estimate on Q s pSq and an estimate on Q f pSq.From the representation of the density in the solid region in Lemma 2.1 (ii) it follows that ρ α Ñ ρ in L 1 pKq for compact sets K Ă Q s pSq and thus osc γ`1 rρ α Ñ ρs pQ s pSqq " 0. In the fluid region the bound is achieved, under exploitation of the effective viscous flux identity (152), by the same arguments as in the all-fluid case, c.f. [16, Proposition 6.1].Finally, the renormalized continuity equation in the limit is also obtained exactly as in the all-fluid case, c.f. [16, Proposition 7.1]: the idea is to pass to the limit in the renormalized continuity equation ( 26) on the α-level for the choice ζ " T k .Thanks to the boundedness of T k , the regularization technique by DiPerna and Lions (c.f.[37, Lemma 6.9]) can be applied to the limit identity.Letting k Ñ 8 and exploiting the bound (153) of the oscillation defect measure, we then obtain the renormalized continuity equation ( 26), ( 27) also for ρ and u.
Having shown the relations (i)-(iii) we now obtain strong convergence of ρ α .Indeed, similarly as in the corresponding relation (142) in the -limit and under exploitation of the concavity of T k , we see that the left-hand side of the effective viscous flux identity (152) is non-negative.This, in combination with the bound (153) of the oscillation defect measure and a comparison between the renormalized continuity equations on the α-level and in the limit, yields, similarly to the first inequality in (142), for τ P r0, T s and ζpsq :" s lnpsq.As in Section 8.3, this inequality implies pointwise convergence of ρ α in p0, T q ˆΩ and therefore z 8 " ρ γ almost everywhere in Q f pSq.In the remaining terms of the momentum equations (143) we may pass to the limit as in the previous limit passages and obtain We are now in the position to conclude the proof of the main result.

Proof of the main result
The function X in ( 21) is defined in (145).The properties of ρ, u and B in ( 22)-( 24), except for the continuity of ρ in time, are shown in ( 146)-(148).The continuity equation (25) and its renormalization ( 26), (27) are derived in (150) and the relation (iii) in Section 9.3, respectively.In particular, ρ as a renormalized solution to the continuity equation satisfies ρ P Cpr0, T s; L 1 pΩqq, c.f. [15, Proposition 4.3], which concludes the proof of (22).The momentum equation ( 28) and the induction equation ( 29) hold true according to (154) and (151).The initial conditions (30) follow as in the previous limit passages, in particular the initial conditions for ρu and B can be derived as the one for the magnetic induction on the η-level in Section 6.2.In the energy inequality (144) on the α-level we can pass to the limit using the weak lower semicontinuity of norms to infer the energy inequality (33).Finally, the compatibility of u with tS i 0 , X i u m i"1 is shown in (149).This finishes the proof of Theorem 2.1.

Appendix
For the construction and estimation of the density on the ∆t-level we use some classical results for the regularized continuity equation which are summarized in the following lemma.
Lemma 10.1 Let T ą 0, r P p0, 1q and assume Ω Ă R 3 to be a bounded domain of class C 2,r .Let V n denote the Galerkin space from our approximate system in Section 3 and let w P Cpr0, T s; V n q.

Proof
The existence of a unique solution ρ to the Neumann problem (155) in the class ρ P C `r0, T s; W 1,2 pΩq ˘č L 2 `0, T ; W 2,2 pΩq ˘, B t ρ P L 2 pp0, T q ˆΩq , which satisfies the estimates (157) and ( 159 l For the limit passage in the time discretization we use a modified version of [38,Theorem 8.9] to infer that different interpolants of the discrete quantities converge to the same limit function: Lemma 10.2 Let Ω Ă R 3 be a bounded domain and assume that f ∆t á f in L 8 p0, T ; L 2 pΩqq, f ∆t á f in L 8 p0, T ; L 2 pΩqq f where f ∆t , f ∆t , f 1 ∆t denote piecewise affine and piecewise constant interpolants of some discrete functions f k ∆t , k " 0, ..., T ∆t as introduced in (59)-(61).Then the limit functions f, f , f 1 coincide, f " f " f 1 a.e. in p0, T q ˆΩ.(160)

Figure 1 :
Figure 1: A domain filled with an electrically conducting compressible fluid and three insulating rigid bodies.
for any b P Y pSq.Finally, for any b P DpΩq with curl b " 0 in a neighbourhood of S 0 , we can argue similarly as in the derivation of the C weak -convergence (113) of B n to deduce that ż Ω B n p¨, xq ¨bpxq dx Ñ ż Ω Bp¨, xq ¨bpxq dx in C pr0, t 0 sq ,

4 3
pQq .The boundary conditions of the limit functions in (123)-(125) follow directly from the corresponding boundary conditions on the η-level, see Proposition 6.1.The set-valued function S in (125) is defined by S :" Sp¨q :" pXp¨; Oqq δ where X, given byX η Ñ X in C `r0, T s; C loc `R3 ˘˘,denotes the solution to the initial value problem dXpt; xq dt " R δ rus pt, Xpt; xqq , Xp0; xq " x @x P R 3 ,

1 βOφ
`θ ρ α are bounded uniformly in L γ`θ pKq and L β`θ pKq, respectively, c.f. [21, Section 4.1], [20, Proposition 2.3].In particular, there exists z 8 P L γ`θ γ pKq such that, after the extraction of a subsequence, In order to identify z 7 , we again need to show strong convergence of ρ α .To this end we use the notion of the oscillation defect measure osc γ`1 rρ α Ñ ρs pOq :" sup kě1 " lim sup αÑ0 ż|T k pρ α q ´Tk pρq| γ`1 dxdt  for measurable sets O Ă p0, T q ˆR3 and a concave cut-off function T k P C 8 pr0, 8qq, k P N, coinciding with the identity function on r0, ks and with 2k on r3k, 8q.The proof of the pointwise convergence of ρ α can be divided into three main steps, each of which consists of showing one of the following three relations, respectively:(i) The effective viscous flux identity ´Tk pρ α q div u α ´Tk pρq div u α ¯dxdt α T k pρ α q ´aρ γ α T k pρq ¯dxdt,
∆t u such that t P rpk ∆t ´1q∆t, k ∆t ∆tq.It holds