Finite Element Procedures for Electroseismic Modeling

Juan E. Santos^{a,b,∗} Fabio I. Zyserman^{a} Patricia M. Gauzellino^{a}

aCONICET, Fac. de Cs. Astronómicas y Geofísicas, Universidad Nacional de La Plata, Paseo del Bosque s/n, B1900FWA - La Plata, Argentina

bDepartment of Mathematics, Purdue University, 150 N. University Street, West Lafayette, Indiana, 47907-2067, USA; santos@math.purdue.edu

∗ Corresponding author, e-mail: santos@math.purdue.edu.ar, fax: 0054-221-4236591 Abstract

This work presents a collection continuous and discrete time finite element procedures for electroseismic modeling in poroelastic fluid-saturated media. The model involves the simultaneous solution of Biot’s equations of motion and Maxwell equations, coupled via an electrokinetc coefficient, employing absorbing boundary conditions at the artificial boundaries. The 2D cases of compressional and vertically polarized shear waves coupled with the transverse magnetic polarization (PSVTMmode) and horizontally polarized shear waves coupled with the transverse electric polarization (SHTE-mode) are first formulated. The finite element procedure for the PSVTM-mode based on rectangular elements is analyzed in detail, including the derivation of a priori error estimates on the continuous-time finite element procedure and stability results for the discrete time procedure. Maxwell equations are discretized in space using the mixed finite element spaces of Nedelec, while while for Biot’ s equations we used a nonconforming element for each component of the solid displacement and the vector part of the Raviart-Thomas-Nedelec of zero order for the fluid displacement vector. The corresponding finite element spaces for the SHTE-mode are later indicated. The finite element spaces and continuous and discrete-time procedures for the full 3D case are formulated and analyzed both for simplicial and parallelepiped elements. Numerical examples are also presented for the 1D SHTE case.

Key words: Seismoelectric Modeling, Poroelasticity, Electromagnetics, Finite element methods,

1 Introduction

cambiar es identico al paper en NMPDE

Seismic waves propagating through near-surface layers of the Earth may induce electromagnetic disturbances that can be measured at the surface (seismoelectric effect)[1–3]. Also, recent tests suggest that the reciprocal process, i.e. surface measurable acoustic disturbances induced by electromagnetic fields (electroseismic effect), is also possible [4,5]. In order to explain these phenomena, Thompson and Gist [6] and Pride [7] suggested that they are generated by an electrokinetic coupling mechanism which can be shortly explained as follows [8,9]. Within a fluid saturated porous medium there exists a nanometer-scale separation of electric charge in which a bound charge existing on the surface of the solid matrix (normally of negative sign) is balanced by adsorbed positive ions of the surrounding fluid, setting an immobile layer. Further from the surface there exists a distribution of mobile counter ions, forming the so called diffuse layer. The effective thickness of this double layer is of about 10 nm. When an electric field is applied to this system, the ions in the diffuse layer move, dragging the pore fluid along with it because of the viscous traction. This is known as electro-osmosis and is responsible for the electroseismic phenomena. On the other hand, the reciprocal situation arises when an applied pressure gradient creates fluid flow and hence, an ionic convection current, which in turn produces an electric field. This is known as electrofiltration and is responsible for the so-called seismoelectric phenomena. Using a volume averaging approach, Pride [7] derived a set of equations describing both electroseismic and seismoelectric effects in electrolyte-saturated porous media. In these equations the coupling mechanism acts through the (generally frequency dependent) electrokinetic coupling coefficient L(ω). When this coefficient is set to zero, Pride’s set of equations turns to the uncoupled Maxwell’s and Biot’s equations, describing the latter mechanical wave propagation in a fluid saturated porous medium [10,11]. There exist already some works implementing different numerical methods to solve the set of equations modeling both mentioned processes. Haartsen and Pride [12], Han and Wang [13], Pain et. al [14], Haines and Pride [9] and White [15] and White and Zhou [16] have proposed several different approaches to numerically study these phenomena.

2 The differential model

Consider a 3D-rectangular domain Ω = Ω_{a}∪Ω_{p} where Ω_{a} and Ω_{p} are associated with the air and subsurface poroelastic (disjoint) parts of Ω, respectively. We will assume that (in cartesian coordinates (x_{1}, x_{2}, x_{3})) all physical quantities describing our domains Ω_{a} and Ω_{p} are independent of the x_{2}-direction (i.e., x_{2} is the symmetry axis) and consider an infinite solenoid J_{m}^{s} in the x_{2}-direction. Including the air region allows us to include in implicit fashion the boundary conditions for the electromagnetic fields at the air-poroelastic interface.

Let us denote by E, H to the electric and magnetic fields in Ω, respectively, and by u^{s}, u^{f} to the solid and relative fluid displacement vectors in Ω_{p}. Under the above symmetry asumption, this source term induces electric and magnetic fields of the form (E_{1}(x_{1}, x_{3}, t), 0, E_{3}(x_{1}, x_{3}, t)), (0, H_{2}(x_{1}, x_{3}, t), 0), respectively, and solid and relative fluid displacements of the form u^{s} = (u_{1}^{s}(x_{1}, x_{3}, t), 0, u_{3}^{s}(x_{1}, x_{3}, t)) and u^{f} = (u_{1}^{f}(x_{1}, x_{3}, t), 0, u_{3}^{f}(x_{1}, x_{3}, t)), respectively. Consequently only compressional and vertically polarized shear seismic waves (PSV-waves) are generated. This is a 2D model known as a PSVTMmode.

Let us identify the 3D vectors (E_{1}(x_{1}, x_{3}, t), 0, E_{3}(x_{1}, x_{3}, t)) and (0, H_{2}(x_{1}, x_{3}, t)), 0) with the 2D vector (E(x_{1}, x_{3}, t) = (E_{1}(x_{1}, x_{3}, t), E_{3}(x_{1}, x_{3}, t)) and the scalar H_{2}(x_{1}, x_{3}, t), respectively. Then recall that

curl H_{2} = (−\frac{∂H_{2}}{∂x_{3}} , \frac{∂H_{2}}{∂x_{1}} ) , curl E = \frac{∂E_{1}}{∂x_{3}} − \frac{∂E_{3}}{∂x_{1}} .

Also, let us identify our 3D-rectangular domain Ω with the 2D-rectangular domain Ω ∩ {y = 0}, so that Ω is the union of the disjoint rectangular subdomains Ω_{a} and Ω_{p}. Let Γ denote the boundary of Ω and let Γ_{a,p} = \overline{Ω}_{a} ∩ \overline{Ω}_{p} denote the free surface. Also let Γ_{a} = ∂Ω_{a} \ Γ_{a,p}, Γ_{p} = ∂Ω_{p} \ Γ_{a,p} denote the artificial boundaries of Ω_{a} and Ω_{p}, respectively.

Following [7,9], for 2D PSVTM electroseismic modeling the electric and magnetic fields E and H and the displacement vectors u^{s} and u^{f} satisfy the coupled electromagnetic-poroelastic equations, stated in the space-time domain as follows:

ε\frac{∂E}{∂t} + σE − curl H_{2} = 0, Ω,  (1)
curl E + \frac{∂H_{2}}{∂t} = J_{m}^{s}, Ω,  (2)
ρ_{b} \frac{∂^{2}u^{s}}{∂t^{2}} + ρ_{f} \frac{∂^{2}u^{f}}{∂t^{2}} − ∇ · τ(u) = 0, Ω_{p},  (3)
ρ_{f} \frac{∂^{2}u^{s}}{∂t^{2}} + m\frac{∂^{2}u^{f}}{∂t^{2}} + \frac{η}{κ_{0}} \frac{∂u^{f}}{∂t} u^{f} − L_{0} \frac{η}{κ_{0}} E + ∇p_{f} = 0, Ω_{p},  (4)
τ_{lm}(u) = 2G ε_{lm}(u^{s}) + δ_{lm} (λ_{c} ∇ · u^{s} + αK_{av} ∇ · u^{f} ) , Ω_{p},  (5)
p_{f}(u) = −αK_{av} ∇ · u^{s} − K_{av}∇ · u^{f}, Ω_{p}.  (6)

In the equations above u = (u^{s}, u^{f}) and τ_{lm}(u) is the stress tensor of the bulk material and p_{f}(u) the fluid pressure, while ε_{lm}(u^{s}) denotes the strain tensor of the solid frame. The coefficients in the stress strain relation are determined as follows. The coefficient G is equal to the elastic shear modulus of the dry matrix. Also,

λ_{c} = K_{c} − G,  (7)

with K_{c} being the bulk modulus of the saturated material. The coefficients in (5)-(6) can be obtained from the relations [34], [35]

α = 1 − \frac{K_{m}}{K_{s}} ,  K_{av} = [\frac{α − ϕ}{K_{s}} + \frac{ϕ}{K_{f}} ]^{−1}  (8)
K_{c} = K_{m} + α^{2}K_{av},

where ϕ denotes the effective porosity and K_{s}, K_{m} and K_{f} denote the bulk modulus of the solid grains composing the solid matrix, the dry matrix and the saturant fluid, respectively.

Also, ε is the electric permitivity, μ the magnetic permeability and σ the conductivity. Furthermore,

ρ_{b} = ϕρ_{f} + (1 − ϕ)ρ_{s},  (9)

where ρ_{s} and ρ_{f} denote the mass densities of the solid grains composing the solid matrix and the saturant fluid. On the other hand, η is the fluid viscosity, κ_{0} the permeability and m is the mass coupling coefficient between the solid and fluid phases in Ω_{p}. Following [?], m can be written in the form

m = \frac{α_{∞}ρ_{f}}{ϕ} ,  (10)

with α_{∞} being the formation tortuosity.

The positive coupling coefficient L_{0} is defined as [12] by

L_{0} = −\frac{ϕ}{α_{∞}} \frac{ϵ_{0}k_{f}ζ}{η} (1 − 2α_{∞} \frac{d̃}{Λ}) ,  (11)

with ζ = 0.008 + 0.026log_{10}(C_{e}) denoting the zeta potential and C_{e} being the electrolyte molarity. In (11) ϵ_{0} and k_{f} are the vacuum and fluid permitivities and

d̃ = \frac{ϵ_{0}k_{f}k_{B}T}{e^{2}z^{2}N_{ic}}  (12)

is the Debye length in meters. in (12) e is the electronic charge, k_{B} is the Boltzman constant, T is the absolute temperature (so that k_{B}T is the thermal energy) z is the ionic valence and N_{ic} the ionic concentration in ions per meters cubed.

To solve equations (1)-(6) in our 2D domain Ω we need a collection of boundary conditions. Let Γ denote the boundary of Ω and let Γ_{a,p} = \overline{Ω}_{a} ∩ \overline{Ω}_{p} denote the free surface. Also let Γ_{a} = ∂Ω_{a} \ Γ_{a,p}, Γ_{p} = ∂Ω_{p} \ Γ_{a,p} denote the artificial boundaries of Ω_{a} and Ω_{p}, respectively. Also, if Γ_{s} is either an inner interface in Ω or a part of the boundaries Γ, Γ_{p} or Γ_{a,p}, set

G_{Γ_{s}}(u) = (τ(u)ν · ν, τ(u)ν · χ, p_{f}(u))^{t},  (13a)
S_{Γ_{s}}(u) = (u^{s} · ν, u^{s} · χ, u^{f} · ν)^{t} ,  (13b)

where , where ^{t} denotes the transpose, ν is the unit outer normal on Γ_{s} and χ is a unit tangent on Γ_{s} oriented counterclockwise.

Then, for ω > 0 consider the solution of (1)-(6) with the absorbing boundary conditions [43], [40]

−ε^{1/2}E · χ + H_{2} = 0, on Γ,  (14)
−G_{Γ_{p}}(u) = DS_{Γ_{p}}(\frac{∂u}{∂t}), on Γ_{p},  (15)

the free surface condition

−G_{Γ_{p}}(u) = 0, on Γ_{a,p},  (16)

and the initial conditions

E(x_{1}, x_{3}, t = 0) = E_{0}, H_{2}(x_{1}, x_{3}, t = 0) = H_{2,0}  (17)

u^{s}(x_{1}, x_{3}, t = 0) = u_{0}^{s}, u^{f}(x_{1}, x_{3}, t = 0) = u_{0}^{f},

\frac{∂u^{s}}{∂t} (x_{1}, x_{3}, t = 0) = u_{1}^{s}, \frac{∂u^{f}}{∂t} (x_{1}, x_{3}, t = 0) = u_{1}^{f}.

The matrix D in (15) is defined as: D = R^{\frac{1}{2}} S^{\frac{1}{2}} R^{\frac{1}{2}} , where S = R^{− \frac{1}{2}} M^{\frac{1}{2}} R^{− \frac{1}{2}} and

R = (\begin{matrix}ρ_{b} & 0 & ρ_{f} \\ 0 & b & 0 \\ ρ_{f} & 0 & ζ \end{matrix}) , M = (\begin{matrix}λ_{c} + 2G & 0 & α K_{av} \\ 0 & G & 0 \\ α K_{av} & 0 & K_{av} \end{matrix}) ,  (18)

where

b = ρ_{b} − \frac{(ρ_{f})^{2}}{m} .

Remark: Note that since α_{∞} ≥ 1, the matrix R is positive definite. Also, we will requiere that the following conditions be satisfied by the coefficients defining the matrix M:

G > 0,  (19a)
λ_{c} + G − α^{2}K_{av} > 0,  (19b)
K_{av} > 0.  (19c)

Conditions (19a), (19b) and (19c) are necessary and sufficient conditions for the matrix M to be positive definite. In particular, the condition (19b) imposes that the inverse of the jacketed compressibility coefficient be strictly positive, see [37]. As a consequence of the positive definitess of the matrices R and M, the matrix D is also positive definite.

3 A weak formulation

For X ⊂ ℝ^{d}, d = 1, 2, 3 with boundary ∂X, let (·, ·)_{X} denote the complex L^{2}(X) inner product for scalar, vector, or matrix valued functions. Also, for s ∈ ℝ, ∥ · ∥_{s,X} will denote the usual norm for the Sobolev space H^{s}(X). In addition, if X = Ω or X = Γ, the subscript X may be omitted such that (·, ·) = (·, ·)_{Ω} or ⟨·, ·⟩ = ⟨·, ·⟩_{Γ}. Set

H(curl, Ω) = {ψ ∈ (L^{2}(Ω))^{2} : curlψ ∈ L^{2}(Ω)},
H(div, Ω_{p}) = {ψ ∈ (L^{2}(Ω_{p}))^{2} : ∇ · ψ ∈ L^{2}(Ω_{p})},
H^{1}(div, Ω_{p}) = {ψ ∈ (H^{1}(Ω_{p}))^{2} : ∇ · ψ ∈ H^{1}(Ω_{p})},

provided with the natural norms

∥ψ∥_{H(curl,Ω)} = (∥ψ∥_{0}^{2} + ∥curlψ∥_{0}^{2})^{\frac{1}{2}} ,
∥ψ∥_{H(div,Ω_{p})} = (∥ψ∥_{0,Ω_{p}}^{2} + ∥∇ · ψ∥_{0,Ω_{p}}^{2})^{\frac{1}{2}} ,
∥ψ∥_{H^{1}(div,Ω_{p})} = (∥ψ∥_{1,Ω_{p}}^{2} + ∥∇ · ψ∥_{1,Ω_{p}}^{2})^{\frac{1}{2}} .

Recall the integration by parts formulas [41,42]

(ψ, curlφ) − (curl ψ, φ) = ⟨ψ · χ, φ⟩ , ∀ψ ∈ H(curl, Ω), φ ∈ H^{1}(Ω),

(20)

(∇ · ψ, φ) + (ψ, ∇φ) = ⟨ψ · ν, φ⟩ , ∀ψ ∈ H(div, Ω), φ ∈ H^{1}(Ω). (21)

To obtain a variational formulation, test (2) against φ ∈ L^{2}(Ω) and test (1) against ψ ∈ H(curl, Ω) and use the integration by parts formula (20) and the boundary condition (14) to obtain the equations

(ε\frac{∂E}{∂t} , ψ) + (σE, ψ) − (H_{2}, curl ψ) + ⟨(\frac{ε}{μ})^{1/2} E · χ, ψ · χ⟩ = 0, (22)
ψ ∈ H(curl, Ω),
(curl E, φ) + (μ\frac{∂H_{2}}{∂t} , φ) = (J_{m}^{s}, φ), φ ∈ L^{2}(Ω).  (23)

Next, test (3) against v^{s} ∈ [H^{1}(Ω_{p})]^{2} and (4) against v^{f} ∈ H(div, Ω_{p}) and use the integration by parts formula (21) and the boundary condition (15). Setting

P = (\begin{matrix}ρ_{b}I_{d} & ρ_{f}I_{d} \\ ρ_{f}I_{d} & mI_{d} \end{matrix})  (24)

where I_{d} is the identity matrix in R^{d} (d = 2 here), we get the equation

(P \frac{∂^{2}u}{∂t^{2}}, v)_{Ω_{p}} + (\frac{η}{κ_{0}} \frac{∂u^{f}}{∂t} , v^{f})_{Ω_{p}} + A(u, v) − (L_{0} \frac{η}{κ_{0}} E, v^{f} )_{Ω_{p}}  (25)

+ ⟨DS_{Γ_{p}}(\frac{∂u}{∂t}), S_{Γ_{p}}(v)⟩_{Γ_{p}} = 0, v = (v^{s}, v^{f}) ∈ [H^{1}(Ω_{p})]^{2} × H(div, Ω_{p}).

In (25) A(u, v) is the bilinear form defined as

A(u, v) = ∑_{l,m} (τ_{lm}(u), ε_{lm}(v^{s}))_{Ω_{p}} − (p_{f}(u), ∇ · v^{f} )_{Ω_{p}} = (M ϵ̃(u), ϵ̃(v))_{Ω_{p}} ,

u, v ∈ [H^{1}(Ω_{p})]^{2} × H(div, Ω_{p}),  (26)

where the matrix M in (26) is given by

M = (\begin{matrix}λ_{c} + 2G & λ_{c} & α K_{av} & 0 \\ λ_{c} & λ_{c} + 2G & α K_{av} & 0 \\ α K_{av} & α K_{av} & K_{av} & 0 \\ 0 & 0 & 0 & 4G \end{matrix}) ,  (27)

and ε̃(u) = (ε_{11}(u^{(s)}), ε_{22}(u^{(s)}), ∇ · u^{(f)}, ε_{12}(u^{(s)})^{t}. Furthermore, we assume that P and M are positive definite since they are associated with the kinetic and strain energy densities, respectively. We also assume that the entries in these two matrices are bounded below and above by positive constants.

Our continuous weak formulation is stated as follows: find (E, H_{2}, u^{s}, u^{f}) ∈ H(curl, Ω) × L^{2}(Ω) × [H^{1}(Ω_{p})]^{2} × H(div, Ω_{p}) satisfying (22), (23) and (25).

Let us analyze the uniqueness of the solution of our continous problem. Set J_{m}^{s} = 0 in (25) and set to zero the initial conditions in (17). Then choose φ = H_{2} in (23) to get

(curl E, H_{2}) + (μ\frac{∂H_{2}}{∂t} , H_{2}) = 0.  (28)

Also, choose ψ = E in (22) and use (28) to obtain

(ε\frac{∂E}{∂t} , E) + (σE, E) + (μ\frac{∂H_{2}}{∂t} , H_{2}) + ⟨(\frac{ε}{μ})^{1/2} E · χ, E · χ⟩ = 0. (29)

Next, choose v^{s} = \frac{∂u^{s}}{∂t} , v^{f} = \frac{∂u^{f}}{∂t} in (25) and add the resulting equation to (29) to get

\frac{1}{2} \frac{d}{dt} [(P \frac{∂u}{∂t}, \frac{∂u}{∂t})_{Ω_{p}} + (M ϵ̃(u), ϵ̃(u))_{Ω_{p}} + (εE, E) + (μH_{2}, H_{2})] (30)

+Φ(E, \frac{∂u^{f}}{∂t}) + (σE, E)_{Ω_{a}} + ⟨(\frac{ε}{μ})^{1/2} E · χ, E · χ⟩

+ ⟨DS_{Γ_{p}}(\frac{∂u}{∂t}), S_{Γ_{p}}(\frac{∂u}{∂t})⟩_{Γ_{p}} = 0,

where

Φ(E, \frac{∂u^{f}}{∂t}) = (σE, E)_{Ω_{p}} − (L_{0} \frac{η}{κ_{0}} E, \frac{∂u^{f}}{∂t} )_{Ω_{p}} + (\frac{η}{κ_{0}} \frac{∂u^{f}}{∂t} , \frac{∂u^{f}}{∂t})_{Ω_{p}}.  (31)

Set

A_{min,p} = inf {A(x_{1}, x_{3}), (x_{1}, x_{3}) ∈ Ω_{p}}, A = σ, κ_{0}, κ_{0,max} = sup {κ_{0}(x_{1}, x_{3}), (x_{1}, x_{3}) ∈ Ω_{p}}.

Assume that

σ_{min,p} > 0, κ_{0,min,p} > 0, κ_{0,max,p} < ∞,  (32)
C_{1} = min (σ_{min,p} − \frac{L_{0}η}{2κ_{0,max}, p}, \frac{η}{κ_{0,max}, p} − \frac{L_{0}η}{2κ_{0,max}, p}) > 0.  (33)

Thus,

Φ(E, \frac{∂u^{f}}{∂t}) ≥ C_{1} (∥E∥_{0,Ω_{p}}^{2} + ∥u^{f}∥_{0,Ω_{p}}^{2}) .  (34)

Next, since the matrix M is positive definite, Korn’ second inequality [51,52] implies that

(M ϵ̃(u), ϵ̃(u))_{Ω_{p}} ≥ C_{2} (∥u^{s}∥_{1,Ω_{p}}^{2} + ∥u^{f}∥_{H(div,Ω_{p}}^{2})  (35)
−C_{3} (∥u^{s}∥_{0,Ω_{p}}^{2} + ∥u^{f}∥_{0,Ω_{p}}^{2}) .

Set

Z = [H^{1}(Ω_{p})]^{2} × H(div, Ω_{p}),

provided with the natural norm

∥u∥_{Z} = (∥u^{s}∥_{1,Ω_{p}}^{2} + ∥u^{f}∥_{H(div,Ω_{p}}^{2}) .

Then choose a constant ζ such that ζ > C_{3} and define the bilinear form

A_{ζ}(u, v) = A(u, v) + ζ(u, v),

so that A_{ζ} is Z-coercive, i.e.,

A_{ζ}(u, u) ≥ C_{2}∥u∥_{Z}^{2},  (36)

Next, add to (30) the inequality

\frac{ζ}{2} \frac{d}{dt}∥u∥_{0,Ω_{p}}^{2} ≤ \frac{ζ}{2} (∥u∥_{0,Ω_{p}}^{2} + ∥\frac{∂u}{∂t}∥_{0,Ω_{p}}^{2}) ,

integrate in time the resulting inequality and use (34) and (36) to obtain

\frac{1}{2} [(P \frac{∂u}{∂t}, \frac{∂u}{∂t})_{Ω_{p}} (t) + C_{2}∥u(t)∥_{Z}^{2} + (εE, E)(t) + (μH_{2}, H_{2})(t)]  (37)

+C_{1} ∫_{0}^{t} (∥E(s)∥_{0,Ω_{p}}^{2} + ∥u^{f}(s)∥_{0,Ω_{p}}^{2}) ds + ∫_{0}^{t}(σE, E)_{Ω_{a}}(s)ds

+ ∫_{0}^{t} ⟨(\frac{ε}{μ})^{1/2} E · χ, E · χ⟩ (s)ds + ∫_{0}^{t} ⟨DS_{Γ_{p}}(\frac{∂u}{∂t}), S_{Γ_{p}}(\frac{∂u}{∂t})⟩_{Γ_{p}} (s)ds

≤ ∫_{0}^{t} (∥u(s)∥_{0,Ω_{p}}^{2} + ∥\frac{∂u}{∂t}(s)∥_{0,Ω_{p}}^{2}) ds.

Assume that ε and μ are bounded below and above by positive constants, so that

0 < ε_{∗} ≤ ε ≤ ε^{∗} < ∞, 0 < μ_{∗} ≤ μ ≤ μ^{∗} < ∞.  (38)

Thus apply Gronwall’s lemma in (37), note that all terms in the left-hand side of (37) are nonnegative and use (38) to conclude that

∥u(t)∥_{Z} = 0, ∥E(t)∥_{0} = 0, ∥H_{2}(t)∥_{0} = 0, ∀t,

so that uniqueness holds for the solution of our differential problem. Assuming sufficient regularity on the initial data and the external sources, existence can be derived using the compacteness argument of Lions [25] with an argument similar to that given in [24]. For brevity the argument is ommited. The result is summarized in the following theorem.

Theorem 1 Assume the validity of (33) and (38) and that the matrices M and D are positive definite. Then there exists a unique solution of problem (2)-(4) with the boundary conditions (14)-(15) and the initial conditions (17).

4 A continuous-time finite element method for the PSVTM-mode. Rectangular case.

Let T^{h}(Ω) be a nonoverlapping quasiregular partition of Ω = Ω_{p} ∪ Ω_{a} into rectangles Ω_{j} of diameter bounded by h such that \overline{Ω} = ∪_{j=1}^{J}\overline{Ω}_{j}. Denote by ξ_{j} and ξ_{jk} the midpoints of Γ_{j} = ∂Ω_{j} ∩Γ and Γ_{jk} = Γ_{kj} = ∂Ω_{j} ∩∂Ω_{k}, respectively.

To approximate the electromagnetic fields E, H_{2} we will employ the mixed finite element space V^{h} × W^{h}, defined as follows [18,22]:

V^{h} = {ψ ∈ H(curl, Ω) : ψ|Ω_{j} ∈ V_{j}^{h} ≡ P_{0,1}(Ω_{j}) × P_{1,0}(Ω_{j})},
W^{h} = {φ ∈ L^{2}(Ω) : φ|Ω_{j} ∈ W_{j}^{h} ≡ P_{0}(Ω_{j})}.

Here, P_{s,t}(Ω_{j}) denote the polynomials of degree not greater than s in x and not greater than t in z on Ω_{j}, while P_{0} denote the constants on Ω_{j}. The functions in V^{h} have continuous tangential components across the internal boundaries Γ_{jk}. Also, curlV^{h} ⊂ W^{(h)} .

Following [22], the degrees of freedom for V^{h} are defined in the following way. Let Ω_{j} be a general element of the partition T^{h}(Ω) and let ψ ∈ [H^{1}(Ω_{j})]^{2}. Then define the following moments on Γ_{jk}:

M_{Γ_{jk}}(ψ) = {⟨ψ · τ, f⟩_{Γ_{jk}} ; f ∈ P_{0}(Γ_{jk})} .  (39)

Note that (39) are curl–conforming and unisolvent for elements in V^{h}.

To approximate each component of the solid displacement vector we employ the nonconforming finite element space N C^{h} as in [23], while to approximate the fluid displacement vector we choose M^{h}, the vector part of the Raviart-Thomas-Nedelec space [17,18] of zero order. More precisely, set

R̂ = [−1, 1]^{2},  N̂ C(R̂) = Span\{1, x̂_{1}, x̂_{3}, α̃(x̂_{1})−α̃(x̂_{3})\}, α̃(x̂_{1}) = x̂_{1}^{2}−\frac{5}{3}x̂_{1}^{4}.

with the degrees of freedom being the values at the midpoint of each edge of R̂. Next, for each Ω_{j} ⊂ Ω_{p}, let F_{Ω_{j}} : R̂ → Ω_{j} be an invertible affine mapping such that F_{Ω_{j}}(R̂) = Ω_{j}, and define

N C_{j}^{h} = {v : v = v̂ ∘ F_{Ω_{j}}^{−1}, v̂ ∈ N̂ C(R̂)}.

Thus,

N C^{h} = {v : v_{j} = v|_{Ω_{j}} ∈ N C_{j}^{h}, v_{j}(ξ_{jk}) = v_{k}(ξ_{jk}) ∀(j, k)},
M^{h} = {w ∈ H(div, Ω_{p}) : w|Ω_{j} ∈ M_{j}^{h} ≡ P_{1,0}(Ω_{j}) × P_{0,1}(Ω_{j})}.

To state the approximating properties of the finite element spaces defined above we introduce the following four projection operators.

First, let

[H_{h}^{1}(Ω)]^{2} = {ψ : ψ|Ω_{j} ∈ [H^{1}(Ω_{j})]^{2}},

with [H_{h}^{1}(Ω_{p})]^{2} defined in similar fashion and if Γ_{jk,p} denotes any inner interface Γ_{jk} in Ω_{p} let

Λ̃^{h} = {λ̃^{h} : λ̃_{jk}^{h} = tr _{Γ_{jk,p}}(λ̃^{h}|_{Ω_{j}}) ∈ [P_{0}(Γ_{jk,p})]^{2} ≡ Λ̃_{jk}^{h}, λ̃_{jk}^{h} + λ̃_{kj}^{h} = 0} ,

where P_{0}(Γ_{jk,p}) denotes the constant functions defined on Γ_{jk,p}.

Remark: Note that there are two copies of [P_{0}(Γ_{jk,p})]^{2} assigned to each Γ_{jk,p}, one from Ω_{j} to Ω_{k} and another from Ω_{k} to Ω_{j}.

Then we define the projections Π_{h} : H(curl, Ω) ∩ [H_{h}^{1}(Ω)]^{2} → V^{h} : ⟨(ψ − Π_{h}ψ) · χ, 1⟩_{B} = 0, B = Γ_{jk} or Γ_{j},

(40)

P_{h} : L^{2}(Ω) → W^{h} : (P_{h}w − w, φ) = 0, φ ∈ W^{h}, (41)
R_{h} : [H^{2}(Ω_{p})]^{2} → [N C^{h}]^{2} :  (v_{i}^{s} − R_{h}v_{i}^{s})(ξ) = 0, ξ = ξ_{jk} or ξ_{j}, (42)
for v^{s} = (v_{1}^{s}, v_{2}^{s}),
Q_{h} : [H^{1}(Ω_{p})]^{2} → M^{h} :  ⟨(v^{f} − Q_{h}v^{f}) · ν, 1⟩_{B} = 0; B = Γ_{jk,p} or Γ_{j}, (43)
S_{h} : [H^{2}(Ω_{p})]^{2} × H^{1}(div; Ω_{p}) → Λ̃^{h} :  ⟨τ(v)ν − S_{h}(v), 1⟩_{B} = 0,  (44)
v = (v^{s}, v^{f}), B = Γ_{jk,p} or Γ_{j}.

Let us define the broken norms

∥v∥_{s,h,Ω_{p}}^{2} = ∑_{Ω_{j}⊂Ω_{p}} ∥v∥_{s,Ω_{j}}^{2}.

The approximation properties of these operators can be stated as follows [18,26]:

∥ψ − Π_{h}ψ∥_{0} ≤ Ch∥ψ∥_{1}, ψ ∈ [H^{1}(Ω)]^{2},  (45)
∥curl(ψ − Π_{h}ψ)∥_{0} ≤ Ch∥curl ψ∥_{1}, ψ ∈ [H^{1}(Ω)]^{2}, curlψ ∈ H^{1}(Ω),  (46)
∥P_{h}φ − φ∥_{0} ≤ Ch∥φ∥_{1},  ∀φ ∈ H^{1}(Ω),  (47)

[∥v^{s} − R_{h}v^{s}∥_{Ω_{p}} + h∥v^{s} − R_{h}v^{s}∥_{1,h,Ω_{p}}^{2} + h^{2}∥v^{s} − R_{h}v^{s}∥_{2,h,Ω_{p}}^{2}

+h^{\frac{1}{2}} ( ∑_{Ω_{j}⊂Ω_{p}} ∥v^{s} − R_{h}v^{s}∥_{0,∂Ω_{j}}^{2} )^{\frac{1}{2}} + h^{\frac{3}{2}} ( ∑_{Ω_{j}⊂Ω_{p}} ∥τ(v_{j})ν_{j} − S_{h}v_{j}∥_{0,∂Ω_{j}}^{2} )^{1/2} (48)

≤ Ch^{2} (∥v^{s}∥_{2,Ω_{p}} + ∥∇ · v^{f}∥_{1,Ω_{p}}) , v = (v^{s}, v^{f} ) ∈ [H^{2}(Ω_{p})]^{2} × H^{1}(div, Ω_{p}),

∥Q_{h}v^{f} − v^{f}∥_{0,Ω_{p}} ≤ Ch∥v^{f}∥_{1,Ω_{p}}, v^{f} ∈ [H^{1}(Ω_{p})]^{2},  (49)
∥∇ · (v^{f} − Q_{h}v^{f})∥_{0,Ω_{p}} ≤ Ch∥∇ · v^{f}∥_{1,Ω_{p}}, v^{f} ∈ H^{1}(div, Ω_{p}).  (50)

Note that since curl ψ ∈ W^{h} ∀ψ ∈ V^{h}, it follows from (41) that

(P_{h}f − f, curl ψ) = 0, ∀ψ ∈ V^{h}.  (51)

Also note the orthogonality property for functions on N C^{h}:

⟨v_{j}^{s} − v_{k}^{s}, 1⟩_{Γ_{jk}} = 0 for all interior interfaces Γ_{jk}, v^{s} ∈ N C^{h}.  (52)

Set

A_{h}(u, v) = ∑_{Ω_{j}⊂Ω_{p}} [∑_{l,m} (τ_{lm}(u), ε_{lm}(v^{(s)}))_{Ω_{j}} − (p_{f}(u), ∇ · v^{f}))_{Ω_{j}} ]  (53)

= ∑_{Ω_{j}⊂Ω_{p}} (M ϵ̃(u), ϵ̃(v))_{Ω_{j}} ,

and

Θ_{h} ((E, H_{2}, u^{s}, u^{f}), (ψ, φ, v^{s}, v^{f})) = (ε\frac{∂E}{∂t} , ψ) + (σE^{,}ψ) − (H_{2}, curl ψ)  (54)

+(curl E, φ) + ⟨(\frac{ε}{μ})^{1/2} E · χ, ψ · χ⟩ + (μ\frac{∂H_{2}}{∂t} , φ)

+ (P \frac{∂^{2}u}{∂t^{2}}, v)_{Ω_{p}} + (\frac{η}{κ_{0}} \frac{∂u^{f}}{∂t} , v^{f} )_{Ω_{p}} + A_{h}(u, v) − (L_{0} \frac{η}{κ_{0}} E, v^{f} )_{Ω_{p}} + ⟨DS_{Γ}(\frac{∂u}{∂t}), S_{Γ}(v)⟩_{Γ_{p}} .

Let J = (0, T) and

Y^{h} = V^{h} × W^{h} × (N C^{h})^{2} × M^{h}.

Then the continuous-time Galerkin procedure is defined as follows: find (E^{h}, H_{2}^{h}, u^{s,h}, u^{f,h}) : J → Y^{h} such that

Θ_{h} ((E^{h}, H_{2}^{h}, u^{s,h}, u^{f,h}), (ψ, φ, v^{s}, v^{f})) = (J_{m}^{s}, φ)  (55)
(ψ, φ, v^{s}, v^{f}) ∈ Y^{h},

with the initial conditions

E^{h}(t = 0) ≈ E_{0}, H_{2}^{h}(t = 0) ≈ H_{2,0}  (56)
u^{s,h}(t = 0) ≈ u_{0}^{s}, u^{f,h}(t = 0) ≈ u_{0}^{f}.  (57)

To analyze the uniqueness for (55)-(57), let us introduce the space

Z^{h} = [H_{h}^{1}(Ω_{p})]^{2} × H(div, Ω_{p})

provided with the norm

∥u∥_{Z^{h}}^{2} = (∥u^{s}∥_{1,h,Ω_{p}}^{2} + ∥u^{f}∥_{H(div,Ω_{p})}^{2})^{1/2} .

As in the continuous case, let us define the bilinear form

A_{ζ,h}(u, v) = A_{h}(u, v) + ζ(u, v), u, v ∈ Z^{h},  (58)

with ζ > C_{3}, so that

A_{ζ,h}(u, u) ≥ C_{2}∥u∥_{Z^{h}}^{2},  (59)
A_{ζ,h}(u, v) ≤ C∥u∥_{Z^{h}}∥v∥_{Z^{h}}.  (60)

Then if u^{h} = (u^{s,h}, u^{f,h}), a repetition of the argument leading to (37) implies that

∥u^{h}∥_{Z^{h}} = 0, ∥E^{h}(t)∥_{0} = 0, ∥H_{2}^{h}(t)∥_{0} = 0, ∀t,

so that uniqueness holds for (55)-(57). Existence follows from finite-dimensionality. Thus we can state the following theorem.

Theorem 2 Under the hypothesis of the positive definitess of the matrices Mand D and the validity of (33) and (38), there exists a unique solution of (55) for every choice of the initial conditions (56)and (57).

5 Apriori error estimates for the PSVTM-mode

Now we derive the error estimate for the procedure (55) stated in the following theorem.

Theorem 3 Assume that the matrices M, D are positive definite and the validity of (33) and (38). Also assume that ε and μ are piecewise constant on Γ. Then the following a priori error estimate holds:

∥E − E^{h}∥_{L^{∞}(J,L^{2}(Ω))} + ∥H_{2} − H_{2}^{h}∥_{L^{∞}(J,L^{2}(Ω))} + ∥u^{s} − u^{s,h}∥_{L^{∞}(J,H_{h}^{1}(Ω_{p}))}

+∥u^{f} − u^{f,h}∥_{L^{∞}(J,H(div,Ω_{p}))} + ∥\frac{∂(u^{s} − u^{s,h})}{∂t})∥_{L^{∞}(J,L^{2}(Ω_{p}))}  (61)

+∥\frac{∂(u^{f} − u^{f,h})}{∂t})∥_{L^{∞}(J,L^{2}(Ω_{p}))} + ∥u^{s} − u^{s,h}∥_{L^{2}(J,L^{2}(Γ_{p}))}

+∥(u^{f} − u^{f,h}) · ν∥_{L^{2}(J,L^{2}(Γ_{p}))}

≤ C [h (N_{0} + N_{1}) + h^{1/2} (M_{0} + M_{1})] ,

where

M_{1}^{2} = ∥u^{s}∥_{L^{∞}(J,H^{2}(Ω_{p}))}^{2} + ∥\frac{∂u^{s}}{∂t} ∥_{L^{∞}(J,H^{2}(Ω_{p}))}^{2} + ∥u^{f}∥_{L^{∞}(J,H^{1}(div,Ω_{p}))}^{2}

+∥\frac{∂u^{f}}{∂t} ∥_{L^{2}(J,H^{1}(div,Ω_{p}))}^{2} + ∥E∥_{L^{2}(J,H^{1}(Ω))}^{2},

N_{1}^{2} = ∥E∥_{L^{∞}(J,H^{1}(Ω))}^{2} + ∥curl E∥_{L^{∞}(J,H^{1}(Ω))}^{2} + ∥\frac{∂E}{∂t} ∥_{L^{∞}(J,H^{1}(Ω))}^{2} + ∥\frac{∂H_{2}}{∂t} ∥_{L^{2}(J,H^{1}(Ω))}^{2},

M_{0} = ∥u^{s}(0)∥_{2,Ω_{p}} + ∥u^{f}(0)∥_{H^{1}(div,Ω_{p})} + ∥\frac{∂u^{s}}{∂t} (0)∥_{2,Ω_{p}} + ∥\frac{∂u^{f}}{∂t} (0)∥_{H^{1}(div,Ω_{p}},

N_{0} = ∥E(0)∥_{1} + ∥H_{2}(0)∥_{1}.

Proof: Set

δ = (E − E^{h}, H_{2} − H_{2}^{h}, u^{s} − u^{s,h}, u^{f} − u^{f,h}) ≡ (δ^{E}, δ_{2}^{H}, δ^{s}, δ^{f} ) ,

γ = (Π_{h}E − E^{h}, P_{h}H_{2} − H_{2}^{h}, R_{h}u^{s} − u^{s,h}, Q_{h}u^{f} − u^{f,h}) ≡ (γ^{E}, γ_{2}^{H}, γ^{s}, γ^{f} ) , and δ^{B} = (δ^{s}, δ^{f}), γ^{B} = (γ^{s}, γ^{f}).

First, use integration by parts element by element and the boundary conditions (14)-(15) to see that, for ψ ∈ H(curl, Ω), φ ∈ L^{2}(Ω) and v = (v^{s}, v^{f}) ∈ [L^{2}(Ω_{p})]^{4} such that v^{s}|_{Ω_{j}} ∈ [H^{1}(Ω_{j})]^{2} and v^{f}|_{Ω_{j}} ∈ H(div, Ω_{j}),

Θ_{h} ((E, H_{2}, u^{s}, u^{f}), (ψ, φ, v^{s}, v^{f})) = (J_{m}^{s}, φ)  (62)

+ ∑_{Ω_{j}⊂Ω_{p}} [⟨τ(u)ν, v^{s}⟩_{∂Ω_{j}\Γ_{p}} − ⟨p_{f}(u), v^{f} · ν⟩_{∂Ω_{j}\Γ_{p}}] .

Now subtract (62) from (55) to see that, for (ψ, φ, v^{s}, v^{f} ) ∈ Y^{h}

Θ_{h} ((δ^{E}, δ_{y}^{H}, δ^{s}, δ^{f}), (ψ, φ, v^{s}, v^{f})) = ∑_{Ω_{j}⊂Ω_{p}} [⟨τ(u)ν, v^{s}⟩_{∂Ω_{j}\Γ_{p}}  (63)
− ⟨p_{f}(u), v^{f} · ν⟩_{∂Ω_{j}\Γ_{p}}] .

Next, by hypothesis u^{s} ∈ [H^{2}(Ω_{p})]^{2}, u^{f} ∈ H^{1}(div, Ω_{p}) and consequently p_{f}(u) ∈ H^{1/2}(∂Ω_{j}). Also, v^{f} ∈ M^{h} ⊂ H(div, Ω) and consequently v^{f} · ν is continuous across the interior interfaces of the elements Ω_{j}. Thus,

∑_{Ω_{j}⊂Ω_{p}} ⟨p_{f}(u), v^{f} · ν⟩_{∂Ω_{j}\Γ_{p}} = 0.  (64)

Also use that v^{s} ∈ N C^{h} is constant on Γ_{jk} and that S_{h}(u) changes sign on each side of Γ_{jk} (c.f. (44)) to see that

∑_{Ω_{j}⊂Ω_{p}} ⟨S_{h}(u), v^{s}⟩_{∂Ω_{j}\Γ_{p}} = 0.  (65)

Thus, (62) becomes

Θ_{h} ((δ^{E}, δ_{y}^{H}, δ^{s}, δ^{f}), (ψ, φ, v^{s}, v^{f})) = ∑_{Ω_{j}⊂Ω_{p}} ⟨τ(u)ν − S_{h}(u), v^{s}⟩_{∂Ω_{j}\Γ_{p}} , (66)
(ψ, φ, v^{s}, v^{f} ) ∈ Y^{h}.

Now (66) yields

Θ_{h} ((γ^{E}, γ_{2}^{H}, γ^{s}, γ^{f}), (ψ, φ, v^{s}, v^{f}))

= Θ_{h} ((Π_{h}E − E, P_{h}H_{2} − H_{2}, R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f} ) , (ψ, φ, v^{s}, v^{f} ))

+ ∑_{Ω_{j}⊂Ω_{p}} ⟨τ(u)ν − S_{h}(u), v^{s}⟩_{∂Ω_{j}\Γ_{p}} , (ψ, φ, v^{s}, v^{f} ) ∈ Y^{h}  (67)

Next, choose ψ = 0, φ = 0, v^{s} = \frac{∂γ^{s}}{∂t} , v^{f} = \frac{∂γ^{f}}{∂t} in (67) and add to the resulting equation the inequality

\frac{ζ}{2} \frac{d}{dt}∥γ^{B}∥_{0,Ω_{p}}^{2} ≤ \frac{ζ}{2} (∥γ^{B}∥_{0,Ω_{p}}^{2} + ∥\frac{∂γ^{B}}{∂t} ∥_{0,Ω_{p}}^{2}) ,

to obtain the inequality

\frac{1}{2} \frac{d}{dt} [(P \frac{∂γ^{B}}{∂t} , \frac{∂γ^{B}}{∂t} )_{Ω_{p}} + A_{h}(γ^{B}, γ^{B})] + (\frac{η}{κ_{0}} \frac{∂γ^{f}}{∂t} , \frac{∂γ^{f}}{∂t} )_{Ω_{p}} + ⟨DS_{Γ_{p}}(\frac{∂γ^{B}}{∂t}), S_{Γ_{p}}(\frac{∂γ^{B}}{∂t})⟩_{Γ_{p}}

≤ \frac{ζ}{2} (∥γ^{B}∥_{0,Ω_{p}}^{2} + ∥\frac{∂γ^{B}}{∂t} ∥_{0,Ω_{p}}^{2})

+ (P(R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f}), \frac{∂γ^{B}}{∂t} )_{Ω_{p}} + A_{h}((R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f}), \frac{∂γ^{B}}{∂t})  (68)

+ (\frac{η}{κ_{0}} \frac{∂(Q_{h}u^{f} − u^{f})}{∂t}, \frac{∂γ^{f}}{∂t} )_{Ω_{p}} + ⟨DS_{Γ_{p}}((R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f})), S_{Γ_{p}}(\frac{∂γ^{B}}{∂t})⟩_{Γ_{p}}

+ (L_{0} \frac{η}{κ_{0}} γ^{E}, \frac{∂γ^{f}}{∂t} )_{Ω_{p}} + (L_{0} \frac{η}{κ_{0}} [Π_{h}E − E], \frac{∂γ^{f}}{∂t} )_{Ω_{p}} + ∑_{Ω_{j}⊂Ω_{p}} ⟨τ(u)ν − S_{h}(u), \frac{∂γ^{s}}{∂t} ⟩_{∂Ω_{j}\Γ_{p}} ,

which will be integrated in time from 0 to t, but first we will bound the time integral of the last seven terms in the right-hand side of (68) First, using the approximation properties (48) and (49)

|(P(R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f}), \frac{∂γ^{B}}{∂t} )_{Ω_{p}} | ≤ C (h^{4}∥u^{s}∥_{2,Ω_{p}}^{2} + h^{2}∥u^{f}∥_{1,Ω_{p}}^{2}

+∥\frac{∂γ^{B}}{∂t} ∥_{2,Ω_{p}}^{2}) ,

so that

∫_{0}^{t} |(P(R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f}), \frac{∂γ^{B}}{∂t} )_{Ω_{p}} (s)| ds  (69)

≤ C [∫_{0}^{t} ∥\frac{∂γ^{B}}{∂t} (s)∥_{0,Ω_{p}}^{2}ds + h^{2}M_{1}] .

Next, using again (49)

∫_{0}^{t} |(\frac{η}{κ_{0}} \frac{∂Q_{h}u^{f} − u^{f})}{∂t}, \frac{∂γ^{f}}{∂t} )_{Ω_{p}} (s)| ds ≤ ϵ̂∫_{0}^{t} (\frac{η}{κ_{0}} \frac{∂γ^{f}}{∂t} , \frac{∂γ^{f}}{∂t} )_{Ω_{p}} (s)ds (70)
+Ch^{2}M_{1}.

Also, using integration by parts in time,

∫_{0}^{t} A_{h}((R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f}), \frac{∂γ^{B}}{∂t})(s)ds  (71)

= A_{h}((R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f}), γ^{B})|_{0}^{t}

− ∫_{0}^{t} A_{h}((\frac{∂R_{h}u^{s} − u^{s}}{∂t} , \frac{∂Q_{h}u^{f} − u^{f}}{∂t}), γ^{B})(s)ds.

Let us bound each term in the right-hand side of (71). First,

|A_{h}((R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f}), γ^{B})(0)|  (72)
≤ C (∥(R_{h}u^{s} − u^{s})(0)∥_{1,h,Ω_{p}} + ∥(Q_{h}u^{f} − u^{f})(0)∥_{H(div,Ω_{p}}) ∥γ^{B}(0)∥_{Z^{h}}
≤ C (∥γ^{B}(0)∥_{Z^{h}}^{2} + h^{2}M_{0}^{2}) .

Now choose the initial condition u^{s,h}(0) ∈ N C^{h}, u^{f,h}(0) ∈ M^{h} such that

A_{ζ,h}(δ^{B}(0), v) = 0, v = (v^{s}, v^{f}) ∈ N C^{h} × M^{h}.  (73)

Then choose v = δ^{B}(0) + (R_{h}u^{s}(0), Q_{h}u^{f}(0)) − (u^{s}(0), u^{f}(0)) in (73) and use the coercivity of A_{ζ,h} in Z^{h} to see that

∥δ^{B}(0)∥_{Z^{h}} ≤ Ch M_{0},  (74)

so that using the triangle inequality and (48) and (49) we obtain the bound

∥γ^{B}(0)∥_{Z^{h}} ≤ Ch M_{0}.  (75)

Thus, (72) becomes

|A_{h}((R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f}), γ^{B})(0)| ≤ Ch^{2}M_{0}^{2}.  (76)

Next,

|A_{h}((R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f}), γ^{B})|(t)|  (77)
≤ ϵ̂∥γ^{B}(t)∥_{Z^{h}}^{2} + Ch^{2} (∥u^{s}(t)∥_{2,Ω_{p}}^{2} + ∥u^{f}(t)∥_{H^{1}(div,Ω_{p})}^{2}) ,

and

∫_{0}^{t} |A_{h}((\frac{∂R_{h}u^{s} − u^{s}}{∂t} , \frac{∂Q_{h}u^{f} − u^{f}}{∂t}), γ^{B})(s)ds|  (78)
≤ C (∫_{0}^{t} ∥γ^{B}(s)∥_{Z^{h}}ds + h^{2}M_{1}^{2}) .

Thus, collecting the estimates (76),(77)and (78) we get the bound

∫_{0}^{t} |A_{h}((R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f}), \frac{∂γ^{B}}{∂t})(s)| ds  (79)
≤ ϵ̂∥γ^{B}(t)∥_{Z^{h}}^{2} + C (∫_{0}^{t} ∥γ^{B}(s)∥_{Z^{h}}ds + h^{2}(M_{0}^{2} + M_{1}^{2})) .

Note that the next term in the right-hand side of (68) can be bounded as follows:

|⟨DS_{Γ_{p}}((R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f})), S_{Γ_{p}}(\frac{∂γ^{B}}{∂t})⟩_{Γ_{p}} |

≤ C ∑_{j} h^{1/2} (∥u^{s}∥_{1/2,Γ_{j}} + ∥u^{f} · ν_{j}∥_{1/2,Γ_{j}} ) (∥\frac{∂γ^{s}}{∂t} ∥_{0,Γ_{j}} + ∥\frac{∂γ^{f} · ν_{j}}{∂t} ∥_{0,Γ_{j}} )

≤ ϵ̂⟨DS_{Γ_{p}}(\frac{∂γ^{B}}{∂t}), S_{Γ_{p}}(\frac{∂γ^{B}}{∂t})⟩_{Γ_{p}} + Ch (∥u^{s}∥_{1,Ω_{p}}^{2} + ∥u^{f}∥_{1,Ω_{p}}^{2}) ,

so that

∫_{0}^{t} |⟨DS_{Γ_{p}}((R_{h}u^{s} − u^{s}, Q_{h}u^{f} − u^{f})), S_{Γ_{p}}(\frac{∂γ^{B}}{∂t})⟩_{Γ_{p}} (s)ds|  (80)

≤ ϵ̂∫_{0}^{t} ⟨DS_{Γ_{p}}(\frac{∂γ^{B}}{∂t}), S_{Γ_{p}}(\frac{∂γ^{B}}{∂t})⟩_{Γ_{p}} (s)ds + Ch M_{1}^{2}.

Next, using (45), the integral in time of the fifth and sixth terms in the righthand side of (68) can be bounded as

∫_{0}^{t} |(L_{0} \frac{η}{κ_{0}} γ^{E}, \frac{∂γ^{f}}{∂t} )_{Ω_{p}} (s)| ds + ∫_{0}^{t} |(L_{0} \frac{η}{κ_{0}} [Π_{h}E − E], \frac{∂γ^{f}}{∂t} )_{Ω_{p}} (s)| ds (81)

≤ ϵ̂∫_{0}^{t} (\frac{η}{κ_{0}} \frac{∂γ^{f}}{∂t} , \frac{∂γ^{f}}{∂t} )_{Ω_{p}} (s)ds + C [h^{2} ∫_{0}^{t} ∥E(s)∥_{1}^{2}ds + ∫_{0}^{t} γ^{E}(s)∥_{0}^{2}ds]

≤ ϵ̂∫_{0}^{t} (\frac{η}{κ_{0}} \frac{∂γ^{f}}{∂t} , \frac{∂γ^{f}}{∂t} )_{Ω_{p}} (s)ds + C [h^{2}M_{1}^{2} + ∫_{0}^{t} γ^{E}(s)∥_{0}^{2}ds]  (82)

Next, using integration by parts in time in the last term in the right-hand side of (68),

∫_{0}^{t} ⟨τ(u)ν − S_{h}(u), \frac{∂γ^{s}}{∂t} ⟩_{∂Ω_{j}\Γ_{p}} (s)ds = ⟨τ(u)ν − S_{h}(u), γ^{s}⟩_{∂Ω_{j}\Γ_{p}} |_{0}^{t}

− ∫_{0}^{t} ⟨τ(\frac{∂u}{∂t})ν − S_{h}(\frac{∂u}{∂t}), γ^{s}⟩_{∂Ω_{j}\Γ_{p}} (s)ds.  (83)

Next, since (τ(u)ν − S_{h}(u))(s) has average value zero on ∂Ω_{j} \ Γ_{p}, if q^{s}(s) is the average value of γ^{s}(s) on Ω_{j} for s = 0, ⋯, t, using the trace inequality and Korn’s second inequality, [51,52],

|⟨τ(u)ν − S_{h}(u), γ^{s}⟩_{∂Ω_{j}\Γ_{p}}( 0)|
= |⟨τ(u)ν − S_{h}(u), γ^{s} − q^{s}⟩_{∂Ω_{j}\Γ_{p}}( 0)|
≤ ∥(τ(u)ν − S_{h}(u))(0)∥_{0,∂Ω_{j}}∥(γ^{s} − q^{s})(0)∥_{0,∂Ω_{j}}  (84)

≤ Ch^{1/2} (∥u^{s}(0)∥_{2,Ω_{j}} + ∥∇ · u^{f}(0)∥_{1,Ω_{j}} ) ∥(γ^{s} − q^{s})(0)∥_{0,Ω_{j}}^{1/2}∥(γ^{s} − q^{s})(0)∥_{1,Ω_{j}}^{1/2}

≤ Ch (∥u^{s}(0)∥_{2,Ω_{j}} + ∥∇ · u^{f}(0)∥_{1,Ω_{j}} ) ∥γ^{s}(0)∥_{1,Ω_{j}}.  (85)

Thus, adding over j in (84) we obtain

∑_{j} |⟨τ(u)ν − S_{h}(u), γ^{s}⟩_{∂Ω_{j}\Γ_{p}} (0)| ≤ Ch^{2}(M_{0}^{2} + M_{1}^{2}).  (86)

Similarly,

|⟨τ(u)ν − S_{h}(u), γ^{s} − q^{s}⟩_{∂Ω_{j}\Γ_{p}} (t)|  (87)

≤ ϵ̂∥γ^{s}(t)∥_{1,Ω_{j}}^{2} + Ch^{2} (∥u^{s}(t)∥_{2,Ω_{j}}^{2} + ∥∇ · u^{f}(t)∥_{1,Ω_{j}}^{2} ) ,

so that adding over j in (87) we get

∑_{j} |⟨τ(u)ν − S_{h}(u), γ^{s} − q^{s}⟩_{∂Ω_{j}\Γ_{p}} (t)|  (88)
≤ ϵ̂∥γ^{s}(t)∥_{1,h}^{2} + Ch^{2}M_{1}^{2}.

In a similar fashion,

∑_{j} |∫_{0}^{t} ⟨τ(\frac{∂u}{∂t})ν − S_{h}(\frac{∂u}{∂t}), γ^{s}⟩_{∂Ω_{j}\Γ_{p}} (s)ds|  (89)
≤ C (∫_{0}^{t} ∥γ^{s}(s)∥_{1,h}^{2}ds + h^{2}M_{1}^{2}) .

Combining the estimates (86), (88) and (89) we conclude that

|∫_{0}^{t} ⟨τ(u)ν − S_{h}(u), \frac{∂γ^{s}}{∂t} ⟩_{∂Ω_{j}\Γ_{p}} (s)ds|  (90)
≤ ϵ̂∥γ^{s}(t)∥_{1,h}^{2} + C (∫_{0}^{t} ∥γ^{s}(s)∥_{1,h}^{2}ds + h^{2}M_{1}^{2}) .

Then add the inequality

\frac{ζ}{2} \frac{d}{dt}∥γ^{B}∥_{0,Ω_{p}}^{2} ≤ \frac{ζ}{2} (∥γ^{B}∥_{0,Ω_{p}}^{2} + ∥\frac{∂γ^{B}}{∂t} ∥_{0,Ω_{p}}^{2}) ,

to (68), integrate the result from 0 to T and use the bounds (69) (70) (79)

(80) and (90) to obtain

\frac{1}{2} (P \frac{∂γ^{B}}{∂t} , \frac{∂γ^{B}}{∂t} )_{Ω_{p}} (t) + A_{ζ,h}(γ^{B}, γ^{B})(t) + ∫_{0}^{t} (\frac{η}{κ_{0}} \frac{∂γ^{f}}{∂t} , \frac{∂γ^{f}}{∂t} )_{Ω_{p}} (s)ds

+ ∫_{0}^{t} ⟨DS_{Γ_{p}}(\frac{∂γ^{B}}{∂t}), S_{Γ_{p}}(\frac{∂γ^{B}}{∂t})⟩_{Γ_{p}} (s)ds  (91)

≤ ϵ̂(∥γ^{B}(t)∥_{Z^{h}}^{2} + ∫_{0}^{t} (\frac{η}{κ_{0}} \frac{∂γ^{f}}{∂t} , \frac{∂γ^{f}}{∂t} )_{Ω_{p}} (s)ds

+ ∫_{0}^{t} ⟨DS_{Γ_{p}}(\frac{∂γ^{B}}{∂t}), S_{Γ_{p}}(\frac{∂γ^{B}}{∂t})⟩_{Γ_{p}} (s)ds)

+C (∫_{0}^{t} (∥γ^{B}(s)∥_{Z^{h}} + ∥\frac{∂γ^{B}}{∂t} (s)∥_{0,Ω_{p}}) ds + h^{2}(M_{0}^{2} + M_{1}^{2}))

+C ∫_{0}^{t} ∥γ^{E}(s)∥_{0}^{2}ds.

Then, use that P and D are positive definite, that A_{ζ,h} is Z^{h}-coercive (see (59)) and apply Gronwall’s lemma in (91) to obtain the following estimate:

∥\frac{∂γ^{B}}{∂t} ∥_{L^{∞}(J,L^{2}(Ω_{p})} + ∥γ^{B}∥_{L^{∞}(J,Z^{h}} + ∥\frac{∂γ^{s}}{∂t} ∥_{L^{2}(J,L^{2}(Γ_{p})}  (92)

+∥\frac{∂γ^{f} · ν}{∂t} ∥_{L^{2}(J,L^{2}(Γ_{p})} ≤ C [h^{1/2}(M_{0} + M_{1}) + ∫_{0}^{t} ∥γ^{E}(s)∥_{0}^{2}ds] .

Next, choose ψ = v^{s} = v^{f} = 0, φ = γ_{2}^{H} in (67) to get

(curl γ^{E}, γ_{2}^{H}) + \frac{1}{2} \frac{d}{dt}(μγ_{2}^{H}, γ_{2}^{H}) = (curl [Π_{h}E − E] , γ_{2}^{H})  (93)

+ (μ\frac{∂(P_{h}H_{2} − H_{2})}{∂t}, γ_{2}^{H}) .

Thus, taking ψ = γ^{E}, φ = 0, v^{s} = v^{f} = 0 in (67) and using (93) yields

\frac{1}{2} \frac{d}{dt} [(εγ^{E}, γ^{E}) + (μγ_{2}^{H}, γ_{2}^{H})] + (σγ^{E}, γ^{E}) + ⟨(\frac{ε}{μ})^{1/2} γ^{E} · χ, γ^{E} · χ⟩

= (\frac{∂(Π_{h}E − E)}{∂t}, γ^{E}) + (curl [Π_{h}E − E] , γ_{2}^{H})  (94)

+ (μ\frac{∂(P_{h}H_{2} − H_{2})}{∂t}, γ_{2}^{H}) + (σ [Π_{h}E − E] , γ^{E}) − ([P_{h}H_{2} − H_{2}] , curl γ^{E})

+ ⟨(\frac{ε}{μ})^{1/2} [Π_{h}E − E] · χ, γ^{E} · χ⟩ .

Let us bound each term in the right-hand side of (94). First note that note that the assumption that the coefficients ε and μ are piecewise contants on Γ, the fact that curl γ^{E} ∈ W^{h} and the orthogonality property (51) that the last two terms in the right-hand side of (94) vanish. The other terms can be bounded using the approximating properties of Π_{h} and P_{h} in (45), (46) and (47) as follows.

|(\frac{∂(Π_{h}E − E)}{∂t}, γ^{E})| + |(curl [Π_{h}E − E] , γ_{2}^{H})|  (95)

+ |(μ\frac{∂(P_{h}H_{2} − H_{2})}{∂t}, γ_{2}^{H})|

≤ C (∥γ_{2}^{E}∥_{0}^{2} + ∥γ_{2}^{H}∥_{0}^{2} + h^{2} (∥curl E∥_{1}^{2} + ∥\frac{∂E}{∂t} ∥_{1}^{2} + ∥\frac{∂H_{2}}{∂t} ∥_{1}^{2})) ..

Thus, apply the bound (95) in (94) to get the inequality

\frac{1}{2} \frac{d}{dt} [(εγ^{E}, γ^{E}) + (μγ_{2}^{H}, γ_{2}^{H})] + (σγ^{E}, γ^{E})  (96)

≤ C (∥γ^{E}∥_{0}^{2} + ∥γ_{2}^{H}∥_{0}^{2} + h^{2} (∥E∥_{1}^{2} + ∥curl E∥_{1}^{2} + ∥\frac{∂E}{∂t} ∥_{1}^{2} + ∥\frac{∂H_{2}}{∂t} ∥_{1}^{2} + ∥\frac{∂u^{f}}{∂t} ∥_{1,Ω_{p}}^{2}))

Now integrate (96) from 0 to t and use (38) to get

∥γ^{E}(t)∥_{0}^{2} + ∥γ_{2}^{H}(t)∥^{2} + ∫_{0}^{t}(σγ^{E}, γ^{E})(s)ds  (97)
≤ C (∥γ^{E}(0)∥_{0}^{2} + ∥γ_{2}^{H}(0)∥_{0}^{2} + ∫_{0}^{t} (γ^{E}(s)∥_{0}^{2} + ∥γ_{2}^{H}(s)∥_{0}^{2}) ds
+Ch^{2}(M_{1}^{2} + N_{1}^{2})) .

Next, choose E^{h}(0) to be the ¶_{h}-projection of E(0) into V^{h} and H_{2}^{h}(0) to be the P_{h}-projection of H_{2}^{h}(0) into W^{h}, so that using the triangle inequality we have that

∥γ^{E}(0)∥ ≤ Ch∥E(0)∥_{1},  ∥γ_{2}^{H}(0)∥ ≤ Ch∥γ_{2}^{H}(0)∥_{1}..  (98)

Thus use (98) and apply Gronwall’s lemma in (97) to obtain the estimate

∥γ^{E}∥_{L^{∞}(J,L^{2}(Ω)} + ∥γ_{2}^{H}∥_{L^{2}(J,L^{2}(Ω)} ≤ Ch(N_{0} + N_{1}).  (99)

Next using (99) in (92) we get

∥\frac{∂γ^{B}}{∂t} ∥_{L^{∞}(J,L^{2}(Ω_{p})} + ∥γ^{B}∥_{L^{∞}(J,Z^{h}} + ∥\frac{∂γ^{s}}{∂t} ∥_{L^{2}(J,L^{2}(Γ_{p})}  (100)

+∥\frac{∂γ^{f} · ν}{∂t} ∥_{L^{2}(J,L^{2}(Γ_{p})} ≤ C [h^{1/2}(M_{0} + M_{1}) + h(N_{0} + N_{1})] .

Finally, using the triangle inequality, the approximating properties (45), (47)

(48) and the estimates (99) and (100) we obtain the estimate in (61). This completes the proof. QED

6 The discrete-time finite element procedures for the PSVTMmode

Let L a positive integer, Δt = T/L, g^{n} = g(nΔt). Set ĝ^{n+1/2} = \frac{g^{n+1}+g^{n}}{2} and

∂^{2}g^{n} = \frac{g^{n+1} − 2g^{n} + g^{n−1}}{Δt^{2}} , ∂g^{n} = \frac{g^{n+1} − g^{n−1}}{2Δt} , d_{t}g^{n} = \frac{g^{n+1} − g^{n}}{Δt} .

Our first fully implicit discrete-time procedure combines a Crank-Nicholson scheme for the Maxwell part and a Galerkin procedure for the Biot-part as follows. Given (E^{h,1}, H_{2}^{h,1}, u^{s,h,0}, u^{s,h,1}, u^{f,h,0}, u^{f,h,1}) ∈ V^{h} × W^{h} × (N C^{h})^{4} ×

(M^{h})^{4}, for n ≥ 1 compute (E^{h,n+1}, H_{2}^{h,n+1}, u^{s,h,n+1}, u^{s,h,n+1}) ∈ Y^{h} such that

(ε d_{t}E^{h,n}, ψ) + (σÊ^{h,n+1/2}, ψ) − (Ĥ_{2}^{h,n+1/2}, curl ψ)  (101)

+ (L_{0} \frac{η}{κ_{0}} ∂u^{f,h,n}, ψ)_{Ω_{p}} + ⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+1/2} · χ, ψ · χ⟩ = 0, ψ ∈ V^{h},

(curlÊ^{h,n+1/2}, φ) + (μd_{t}H_{2}^{h,n}, φ) = 0, φ ∈ W^{h},  (102)
(P∂^{2}u^{h,n}, v)_{Ω_{p}} + (\frac{η}{κ_{0}} ∂u^{f,h,n}, v^{f} )_{Ω_{p}} + A_{h}(\frac{u^{h,n+1} + u^{h,n−1}}{2} , v)  (103)
+ ⟨DS_{Γ}(∂u^{h,n}), S_{Γ}(v)⟩_{Γ_{p}} = (F^{n}, v), v = (v^{s}, v^{f}) ∈ (N C^{h})^{2} × M^{h}.

To analyze the procedure (101)-(103), choose ψ = Ê^{h,n+1/2} in (101) and φ = Ĥ_{2}^{h,n+1/2} in (102) and add the resulting equations to get

\frac{1}{2Δt} [(ε E^{h,n+1}, E^{h,n+1}) − (ε E^{h,n}, E^{h,n}) + (μH_{2}^{h,n+1}, H_{2}^{h,n+1}) − (μH_{2}^{h,n}, H_{2}^{h,n})]
+ (L_{0} \frac{η}{κ_{0}} ∂u^{f,h,n}, E^{h,n+1})_{Ω_{p}} + ⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+1/2} · χ, E^{h,n+1} · χ⟩ = 0. (104)

Next, choose v = ∂u^{h,n} in (103) and add the resulting equation to (104) to obtain

\frac{1}{2Δt} [(ε E^{h,n+1}, E^{h,n+1}) − (ε E^{h,n}, E^{h,n}) + (μH_{2}^{h,n+1}, H_{2}^{h,n+1})  (105)

−(μH_{2}^{h,n}, H_{2}^{h,n}) + (Pd_{t}u^{h,n}, d_{t}u^{h,n})_{Ω_{p}} − (Pd_{t}u^{h,n−1}, d_{t}u^{h,n−1})_{Ω_{p}}

+\frac{1}{4}A_{h}(u^{h,n+1}, u^{h,n+1}) − \frac{1}{4}A_{h}(u^{h,n−1}, u^{h,n−1})]

+Φ(Ê^{h,n+1/2}, ∂u^{f,h,n}) + (σÊ^{h,n+1/2}, Ê^{h,n+1/2})_{Ω_{a}}

+ ⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+1/2} · χ, E^{h,n+1/2} · χ⟩ + ⟨DS_{Γ}(∂u^{h,n}), S_{Γ}(∂u^{h,n})⟩_{Γ_{p}} = (F^{n}, ∂u^{h,n}),

where Φ(Ê^{h,n+1/2}, ∂u^{f,h,n}) is defined in (31). Add the inequality

\frac{ζ}{4Δt} (∥u^{h,n+1}∥_{0.Ω_{p}}^{2} − ∥u^{h,n−1}∥_{0.Ω_{p}}^{2})
≤ C (∥u^{h,n+1}∥_{0.Ω_{p}}^{2} + ∥u^{h,n}∥_{0.Ω_{p}}^{2} + ∥u^{h,n−1}∥_{0.Ω_{p}}^{2}  (106)
+∥d_{t}u^{h,n}∥_{0.Ω_{p}}^{2} + ∥d_{t}u^{h,n−1}∥_{0.Ω_{p}}^{2})

to (105), multiply the result by Δt, add from n = 1 to n = N and use (34) and (59) to obtain

(ε E^{h,N+1}, E^{h,N+1}) + (μH_{2}^{h,N+1}, H_{2}^{h,N+1}) + (Pd_{t}u^{h,N}, d_{t}u^{h,N} )_{Ω_{p}}  (107)

+\frac{C_{2}}{4} (∥u^{h,N+1}∥_{Z^{h}}^{2} + ∥u^{h,N}∥_{Z^{h}}^{2}) + C_{1} ∑_{n=1}^{N} (∥Ê^{h,n+1/2}∥_{0,Ω_{p}}^{2} + ∥∂u^{f,h,n})∥_{0,Ω_{p}}^{2}) Δt

+ ∑_{n=1}^{N} ((σÊ^{h,n+1/2}, Ê^{h,n+1/2})_{Ω_{a}} + ⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+1/2} · χ, E^{h,n+1} · χ⟩

+ ⟨DS_{Γ}(∂u^{h,n}), S_{Γ}(∂u^{h,n})⟩_{Γ_{p}}) Δt

≤ C (∥E^{h,1}∥_{0}^{2} + ∥H_{2}^{h,1}∥_{0}^{2} + ∥u^{h,0}∥_{Z^{h}}^{2} + ∥u^{h,1}∥_{Z^{h}}^{2}

+∥d_{t}u^{h,0}∥_{0,Ω_{p}}^{2} + ∑_{n=1}^{N} (∥F^{n}∥_{0,Ω_{p}}^{2} + +∥u^{h,n+1}∥_{0,Ω_{p}}^{2} + ∥u^{h,n}∥_{0,Ω_{p}}^{2} + ∥u^{h,n−1}∥_{0,Ω_{p}}^{2}

+∥d_{t}u^{h,n}∥_{0,Ω_{p}}^{2} + ∥d_{t}u^{h,n−1}∥_{0,Ω_{p}}^{2}) Δt) .

Next apply Gronwall’s lemma in (107), use (38) and that the matrices P and D are positive definite to conclude that

max_{1≤n≤L−1} [∥E^{h,n}∥_{0}^{2} + ∥H^{h,n}∥_{0}^{2} + ∥d_{t}u^{h,n}∥_{0,Ω_{p}}^{2} + ∥u^{h,n}∥_{Z^{h}}^{2}]  (108)

+ ∑_{n=1}^{L−1} (∥Ê^{h,n+1/2}∥_{0,Ω_{p}}^{2} + ∥∂u^{f,h,n}∥_{0,Ω_{p}}^{2}Δt)

∑_{n=1}^{L−1} ((σÊ^{h,n+1/2}, Ê^{h,n+1/2})_{Ω_{a}} + ∥Ê^{h,n+1/2} · χ∥_{0,Γ}^{2}

+∥∂u^{s,h,n}∥_{0,Γ_{p}}^{2} + ∥∂u^{f,h,n} · ν∥_{0,Γ_{p}}^{2}) Δt

≤ C (∥E^{h,1}∥_{0} + ∥H_{2}^{h,1}∥_{0} + ∥u^{h,0}∥_{Z^{h}} + ∥u^{h,1}∥_{Z^{h}} + ∥d_{t}u^{h,0}∥_{0,Ω_{p}}

+ ∑_{n=1}^{L−1} ∥F^{n}∥_{0,Ω_{p}}^{2}Δt) .  (109)

The estimate (108) yields uniqueness, existence and unconditional stability for the procedure (101)-(103).

A second fully implicit in time procedure combines a Backward Euler in time for Maxwell equations with the implicit procedure used for Biot’s equations in (103). It can be formulated as follows: Given (E^{h,1}, H_{2}^{h,1}, u^{s,h,0}, u^{s,h,1}, u^{f,h,0}, u^{f,h,1}) ∈

V^{h}×W^{h}×(N C^{h})^{4}×(M^{h})^{4}, for n ≥ 1 compute (E^{h,n+1}, H_{2}^{h,n+1}, u^{s,h,n+1}, u^{s,h,n+1}) ∈ Y^{h} such that

(ε d_{t}E^{h,n}, ψ) + (σÊ^{h,n+1}, ψ) − (Ĥ_{2}^{h,n+1/2}, curl ψ)  (110)

+ (L_{0} \frac{η}{κ_{0}} ∂u^{f,h,n}, ψ)_{Ω_{p}} + ⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+1} · χ, ψ · χ⟩ = 0, ψ ∈ V^{h},

(curlÊ^{h,n+1}, φ) + (μd_{t}H_{2}^{h,n}, φ) = 0, φ ∈ W^{h},  (111)
(P∂^{2}u^{h,n}, v)_{Ω_{p}} + (\frac{η}{κ_{0}} ∂u^{f,h,n}, v^{f} )_{Ω_{p}} + A_{h}(\frac{u^{h,n+1} + u^{h,n−1}}{2} , v)  (112)
+ ⟨DS_{Γ}(∂u^{h,n}), S_{Γ}(v)⟩_{Γ_{p}} = (F^{n}, v), v = (v^{s}, v^{f}) ∈ (N C^{h})^{2} × M^{h}.

This procedure can be analyzed with a similar argument to that given for the procedure (101)-(103). The analysis, omitted here for brevity, yields existence, uniqueness and unconditional stability for this method.

Our third discrete-time procedure consists of Leapfrog scheme for the Maxwell part with the same implicit-time procedure used above for the Biot-part and is formulated as follows. Given (E^{h,0}, H_{2}^{h,1/2}, u^{s,h,0}, u^{s,h,1}, u^{f,h,0}, u^{f,h,1}) ∈ V^{h} ×

W^{h} ×(N C^{h})^{4} ×(M^{h})^{4}, for n ≥ 1 compute (E^{h,n}, H_{2}^{h,n+1/2}, u^{s,h,n+1}, u^{s,h,n+1}) ∈ Y^{h} such that

(ε d_{t}E^{h,n}, ψ) + (σÊ^{h,n+1/2}, ψ) − (H_{2}^{h,n+1/2}, curl ψ)  (113)

+ (L_{0} \frac{η}{κ_{0}} ∂u^{f,h,n}, ψ)_{Ω_{p}} + ⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+1/2} · χ, ψ · χ⟩ = 0, ψ ∈ V^{h},

(curlÊ^{h,n+1/2}, φ) + (μ\frac{H_{2}^{h,n+3/2} − H_{2}^{h,n+1/2}}{Δt} , φ) = 0, φ ∈ W^{h},  (114)
(P∂^{2}u^{h,n}, v)_{Ω_{p}} + (\frac{η}{κ_{0}} ∂u^{f,h,n}, v^{f} )_{Ω_{p}} + A_{h}(\frac{u^{h,n+1} + u^{h,n−1}}{2} , v)  (115)
+ ⟨DS_{Γ}(∂u^{h,n}), S_{Γ}(v)⟩_{Γ_{p}} = (F^{n}, v), v = (v^{s}, v^{f}) ∈ (N C^{h})^{2} × M^{h}.

To analyze the procedure (113)-(115), choose ψ = Ê^{h,n+1/2} in (113) and φ = H_{2}^{h,n+1/2} in (114), add the resulting equations to obtain

\frac{1}{2Δt} [(εE^{h,n+1}, E^{h,n+1}) − (εE^{h,n}, E^{h,n})] + (σÊ^{h,n+1/2}, Ê^{h,n+1/2})  (116)

+(μ\frac{H_{2}^{h,n+3/2} − H_{2}^{h,n+1/2}}{Δt} , H_{2}^{h,n+3/2})

−(μ\frac{H_{2}^{h,n+3/2} − H_{2}^{h,n+1/2}}{Δt} , H_{2}^{h,n+3/2} − H_{2}^{h,n+1/2}) + (L_{0} \frac{η}{κ_{0}} ∂u^{f,h,n}, Ê^{h,n+1/2})_{Ω_{p}}

+ ⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+1/2} · χ, E^{h,n+1/2} · χ⟩ = 0.

Now add (116) and the corresponding equation to (116) for the case n = n+1 and use that (114) yields

curl Ê^{h,n+1/2} = −μ (\frac{H_{2}^{h,n+3/2} − H_{2}^{h,n+1/2}}{Δt} )  (117)

Fabio: fijate si las unidades de la ecuacion (117) estan bien, para obtener esa ecuacion dije lo siguiente: tomo funcion de prueba

φ = curl Ê^{h,n+1/2} + μ (\frac{H_{2}^{h,n+3/2} − H_{2}^{h,n+1/2}}{Δt} )

en (114) para obtener

∥curl Ê^{h,n+1/2} + μ (\frac{H_{2}^{h,n+3/2} − H_{2}^{h,n+1/2}}{Δt} ) ∥_{0}^{2} = 0

Como la norma es cero, la funcion es cero en L^{2}(Ω) y entonces

curl Ê^{h,n+1/2} + μ (\frac{H_{2}^{h,n+3/2} − H_{2}^{h,n+1/2}}{Δt} ) = 0

que es lo que dice arriba en (117). No les pongo nro para no confundirte mas

Esta es la manera en que aparece este termino (curl, curl) en el miembro derecho de la (118) debajo, que es el que trae el problema de unidades.

to get the equation

\frac{1}{2Δt} [(εE^{h,n+2}, E^{h,n+2}) − (εE^{h,n}, E^{h,n})] + (σÊ^{h,n+3/2}, Ê^{h,n+3/2}) (118)

+(σÊ^{h,n+1/2}, Ê^{h,n+1/2}) + (μ\frac{H_{2}^{h,n+5/2} − H_{2}^{h,n+1/2}}{Δt} , H_{2}^{h,n+3/2})

+ (L_{0} \frac{η}{κ_{0}} ∂u^{f,h,n}, Ê^{h,n+1/2})_{Ω_{p}} + (L_{0} \frac{η}{κ_{0}} ∂u^{f,h,n+1}, Ê^{h,n+3/2})_{Ω_{p}}

+ ⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+3/2} · χ, E^{h,n+3/2} · χ⟩

+ ⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+1/2} · χ, E^{h,n+1/2} · χ⟩

= (curl E^{h,n+1/2}, \frac{Δt}{μ}curl E^{h,n+1/2}) .

Next add (118) to the two equations obtained by choosing v = ∂u^{h,n} and v = ∂u^{h,n+1} in (115) to get

\frac{1}{2Δt} [(ε E^{h,n+2}, E^{h,n+2}) − (ε E^{h,n}, E^{h,n}) + (Pd_{t}u^{h,n+1}, d_{t}u^{h,n+1})_{Ω_{p}}  (119)

− (Pd_{t}u^{h,n−1}, d_{t}u^{h,n−1})_{Ω_{p}} + \frac{1}{4}A_{h}(u^{h,n+2}, u^{h,n+2}) + \frac{1}{4}A_{h}(u^{h,n+1}, u^{h,n+1})

−\frac{1}{4}A_{h}(u^{h,n}, u^{h,n}) − \frac{1}{4}A_{h}(u^{h,n−1}, u^{h,n−1})] + (μ\frac{H_{2}^{h,n+5/2} − H_{2}^{h,n+1/2}}{Δt} , H_{2}^{h,n+3/2})

+Φ(Ê^{h,n+3/2}, ∂u^{f,h,n+1}) + +Φ(Ê^{h,n+1/2}, ∂u^{f,h,n}) + (σÊ^{h,n+3/2}, Ê^{h,n+3/2})_{Ω_{a}}

+(σÊ^{h,n+1/2}, Ê^{h,n+1/2})_{Ω_{a}} + ⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+3/2} · χ, E^{h,n+3/2} · χ⟩

+ ⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+1/2} · χ, E^{h,n+1} · χ⟩ + ⟨DS_{Γ}(∂u^{h,n+1}), S_{Γ}(∂u^{h,n+1})⟩_{Γ_{p}}

+ ⟨DS_{Γ}(∂u^{h,n}), S_{Γ}(∂u^{h,n})⟩_{Γ_{p}}

= (F^{n+1}, ∂u^{h,n+1}) + (F^{n}, ∂u^{h,n}) + (curl E^{h,n+1/2}, \frac{Δt}{μ}curl E^{h,n+1/2}) .

Next, add (106) and the corresponding inequality for n = n + 1 to (119), multiply the resulting inequality by Δt and sum from n = 1 to n = N and use (34) and (59) to get

\frac{1}{2} [(εE^{h,N+2}, E^{h,N+2}) + (εE^{h,N+1}, E^{h,N+1})] + (Pd_{t}u^{h,N+1}, d_{t}u^{h,N+1})_{Ω_{p}}

+ (Pd_{t}u^{h,N}, d_{t}u^{h,N} )_{Ω_{p}} + \frac{C_{2}}{4} (∥u^{h,N+2}∥_{Z^{h}}^{2} + ∥u^{h,N+1}∥_{Z^{h}}^{2} + ∥u^{h,N}∥_{Z^{h}}^{2})

+ ∑_{n=1}^{N+1} C_{1} (ˆ∥E^{h,n+1/2}∥_{0,Ω_{P}}^{2} + ∥∂u^{f,h,n}∥_{0,Ω_{p}}^{2}) Δt + ∑_{n=1}^{N+1}(Ê^{h,n+1/2}, Ê^{h,n+1/2})_{Ω_{a}}Δt

+ ∑_{n=1}^{N+1} (⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+1/2} · χ, Ê^{h,n+1/2} · χ⟩  (120)

+ ⟨DS_{Γ}(∂u^{h,n}), S_{Γ}(∂u^{h,n})⟩_{Γ_{p}}) Δt + (μH_{2}^{h,N+5/2}, H_{2}^{h,N+3/2})

≤ C_{8} [∥E^{h,2}∥_{0}^{2} + ∥E^{h,1}∥_{0}^{2} + (μH_{2}^{h,3/2}, H_{2}^{h,1/2}) + ∥d_{t}u^{h,0}∥_{0,Ω_{p}}^{2}

+∥d_{t}u^{h,1}∥_{0,Ω_{p}}^{2} + ∥u^{h,0}∥_{Z^{h}}^{2} + ∥u^{h,1}∥_{Z^{h}}^{2} + ∥u^{h,2}∥_{Z^{h}}^{2}

+ ∑_{n=1}^{N} (∥F^{n}∥_{0,Ω_{p}}^{2} + ∥F^{n+1}∥_{0,Ω_{p}}^{2} + (Pd_{t}u^{h,n+1}, d_{t}u^{h,n+1})_{Ω_{p}}

+ (Pd_{t}u^{h,n}, d_{t}u^{h,n})_{Ω_{p}} + (Pd_{t}u^{h,n−1}, d_{t}u^{h,n−1})_{Ω_{p}}

+∥u^{h,n+1}∥_{0,Ω_{p}}^{2} + ∥u^{h,n}∥_{0,Ω_{p}}^{2}) Δt]

+ ∑_{n=1}^{N} Δt (\frac{1}{μ}curl E^{h,n+1/2}, curl E^{h,n+1/2}) Δt.

Fabio: lo que dije en el mail, mas claro aqui: en la desigualdad de arriba, las unidades de

\frac{1}{2} [(εE^{h,N+2}, E^{h,N+2}) + (εE^{h,N+1}, E^{h,N+1})]

deberian ser las mismas que las de

∑_{n=1}^{N} Δt (\frac{1}{μ}curl E^{h,n+1/2}, curl E^{h,n+1/2}) Δt.

Let us bound the last term in the right-hand side of (120). By the quasiregularity assumption on the finite element partition T^{h}, the exists a positive constant C_{4} (C_{4} = \sqrt{4}8) such that

∥curl E^{h,n+1/2}∥_{0} ≤ C_{4}h^{−1}∥E^{h,n+1/2}∥_{0}.  (121)

Thus, using (38),

(\frac{1}{μ}curl E^{h,n+1/2}, curl E^{h,n+1/2})  (122)

≤ \frac{C_{4}^{2}h^{−2}}{2ε^{∗}μ_{∗}} ((εE^{h,n+1}, E^{h,n+1}) + (εE^{h,n+1}, E^{h,n+1})) .

Thus we conclude that

∑_{n=1}^{N} Δt (\frac{1}{μ}curl E^{h,n+1/2}, curl E^{h,n+1/2}) Δt  (123)

≤ ∑_{n=1}^{N} Δt\frac{C_{4}^{2}h^{−2}}{2ε^{∗}μ_{∗}} ((εE^{h,n+1}, E^{h,n+1}) + (εE^{h,n+1}, E^{h,n+1})) Δt

≤ C_{8} ∑_{n=1}^{N} ((εE^{h,n+1}, E^{h,n+1}) + (εE^{h,n}, E^{h,n})) Δt,

provided

Δt\frac{C_{4}^{2}h^{−2}}{2ε^{∗}μ_{∗}} ≤ C_{8}  (124)

or equivalently, assume the stability constraint

Δt ≤ C_{8} \frac{2ε^{∗}μ_{∗}h^{2}}{C_{4}^{2}} .  (125)

Fabio: la desigualdad (125) es de la forma clasica Δt ≤ Ch^{2} valida para metods explicitos. ES posible que las unidades en la (123) esten bien (yo creo que si, al menos uno comienza con la (122) que tiene bien las unidades).

Thus we conclude that

∑_{n=1}^{N} Δt (\frac{1}{μ}curl E^{h,n+1/2}, curl E^{h,n+1/2}) Δt  (126)

≤ C_{8} ∑_{n=1}^{N} ((εE^{h,n+1}, E^{h,n+1}) + (εE^{h,n}, E^{h,n})) Δt.

Next we analyze the last term in the left-hand side of (120). First, right this term in the form

(μH_{2}^{h,N+5/2}, H_{2}^{h,N+3/2}) = (μ [H_{2}^{h,N+5/2} − H_{2}^{h,N+3/2})] , H_{2}^{h,N+3/2}) (127)
+(μH_{2}^{h,N+3/2}, H_{2}^{h,N+3/2}) ≡ T_{1} + T_{2}.

Note that from (114) and (121)

Fabio aca uso lo de siempre, que para el ϵ que yo quiera,

ab ≤ \frac{1}{2}(ϵa^{2} + \frac{1}{ϵ}b^{2})

|T_{1}| = |−Δt (curl E^{h,N+3/2}, H_{2}^{h,N+3/2})| ≤ Δt∥curl E^{h,N+3/2}∥_{0}∥H_{2}^{h,N+3/2}∥_{0}

≤ \frac{μ_{∗}}{2} ∥H_{2}^{h,N+3/2}∥_{0}^{2} + (Δt)^{2} \frac{C_{4}^{2}h^{−2}}{2μ_{∗}} ∥E^{h,N+3/2}∥_{0}^{2}

≤ \frac{μ_{∗}}{2} ∥H_{2}^{h,N+3/2}∥_{0}^{2} + (Δt)^{2} \frac{C_{4}^{2}h^{−2}}{4ε^{∗}μ_{∗}} ((εE^{h,n+1}, E^{h,n+1}) + (εE^{h,n+1}, E^{h,n+1})) .

Assume that

(Δt)^{2} \frac{C_{4}^{2}h^{−2}}{4ε^{∗}μ_{∗}} ≤ \frac{1}{4}

so that, in addition to (125), Δt satifies the constraint

Δt ≤ \frac{(μ_{∗}ε^{∗})^{1/2}}{C_{4}} h.  (128)

Fabio: esto parece estar bien, habia un delta t dividiendo que estaba mal

Lo que estoy haciendo es hacer que este termino sea la mitad del termino

\frac{1}{2} [(εE^{h,N+2}, E^{h,N+2}) + (εE^{h,N+1}, E^{h,N+1})]

por eso pido el \frac{1}{4}. En otras palabras, lo que estoy usando es que

(Δt)^{2} \frac{C_{4}^{2}h^{−2}}{4ε^{∗}μ_{∗}} es adimensional, si no fuera asi no lo puedo restar al \frac{1}{2}.

Si esto esta bien, veamos como arreglar o entender la primera condicion de estabilidad (125), que es practicamente el cuadrado de esta salvo un factor de 2. Donde carajo esta el error en la derivacion de la (125)?

Then from (127) and (128), since

|T_{2}| ≥ μ_{∗}∥H_{2}^{h,N+3/2}∥_{0}^{2},

(μH_{2}^{h,N+5/2}, H_{2}^{h,N+3/2}) ≥ \frac{μ_{∗}}{2} ∥H_{2}^{h,N+3/2}∥_{0}^{2}  (129)

−\frac{1}{4} [((εE^{h,N+2}, E^{h,N+2}) + (εE^{h,N+1}, E^{h,n+1}))] .

Using (126) and (129) in (120) we obtain

\frac{1}{4} [(εE^{h,N+2}, E^{h,N+2}) + (εE^{h,N+1}, E^{h,N+1})] + \frac{μ_{∗}}{2} ∥H_{2}^{h,N+3/2}∥_{0}^{2}  (130)

+ (Pd_{t}u^{h,N+1}, d_{t}u^{h,N+1})_{Ω_{p}} + (Pd_{t}u^{h,N}, d_{t}u^{h,N} )_{Ω_{p}}

+\frac{C_{2}}{4} (∥u^{h,N+2}∥_{Z^{h}}^{2} + ∥u^{h,N+1}∥_{Z^{h}}^{2} + ∥u^{h,N}∥_{Z^{h}}^{2})

+ ∑_{n=1}^{N+1} [C_{1} (Ê^{h,n+1/2}∥_{0,Ω_{P}}^{2} + ∥∂u^{f,h,n}∥_{0,Ω_{p}}) + ∥Ê^{h,n+1/2}∥_{0,Ω_{a}}] Δt (131)

+ ∑_{n=1}^{N+1} (⟨(\frac{ε}{μ})^{1/2} Ê^{h,n+1/2} · χ, Ê^{h,n+1/2} · χ⟩

+ ⟨DS_{Γ}(∂u^{h,n}), S_{Γ}(∂u^{h,n})⟩_{Γ_{p}}) Δt

≤ C [∥E^{h,2}∥_{0}^{2} + ∥E^{h,1}∥_{0}^{2} + ∥H_{2}^{h,1/2}∥_{0}^{2} + ∥H_{2}^{h,3/2}∥_{0}^{2} + ∥d_{t}u^{h,0}∥_{0,Ω_{p}}^{2}

+∥d_{t}u^{h,1}∥_{0,Ω_{p}}^{2} + ∥u^{h,0}∥_{Z^{h}}^{2} + ∥u^{h,1}∥_{Z^{h}}^{2} + ∥u^{h,2}∥_{Z^{h}}^{2}

∑_{n=1}^{N} (∥F^{n}∥_{0,Ω_{p}}^{2} + ∥F^{n+1}∥_{0,Ω_{p}}^{2} + (εE^{h,n+2}, E^{h,n+2})
+(εE^{h,n+1}, E^{h,n+1}) + (Pd_{t}u^{h,n+1}, d_{t}u^{h,n+1})_{Ω_{p}}
+ (Pd_{t}u^{h,n}, d_{t}u^{h,n})_{Ω_{p}} + (Pd_{t}u^{h,n−1}, d_{t}u^{h,n−1})_{Ω_{p}}
+∥u^{h,n+1}∥_{0,Ω_{p}}^{2} + ∥u^{h,n}∥_{0,Ω_{p}}^{2}) Δt] .

Next apply the discrete Gronwall’s lemma in (119), use the positive definitess of the matrices P and D and (38) to derive the estimate

max_{1≤n≤L−1} [∥E^{h,n}∥_{0}^{2} + ∥H_{2}^{h,n+1/2}∥_{0}^{2} + ∥d_{t}u^{h,n}∥_{0,Ω_{p}}^{2} + ∥u^{h,n}∥_{Z^{h}}^{2}]  (132)

+ ∑_{n=1}^{L−1} (∥Ê^{h,n}∥_{0,Ω_{p}}^{2} + ∥∂u^{f,h,n}∥_{0,Ω_{p}}^{2}Δt)

∑_{n=1}^{L−1} ((σÊ^{h,n+1/2}, Ê^{h,n+1/2})_{Ω_{a}} + ∥Ê^{h,n+1/2} · χ∥_{0,Γ}^{2}

+∥∂u^{s,h,n}∥_{0,Γ_{p}}^{2} + ∥∂u^{f,h,n} · ν∥_{0,Γ_{p}}^{2}) Δt

≤ C (∥E^{h,1}∥_{0} + ∥E^{h,2}∥_{0} + ∥H_{2}^{h,1/2}∥_{0} + ∥H_{2}^{h,3/2}∥_{0} + ∥u^{h,0}∥_{Z^{h}} + ∥u^{h,1}∥_{Z^{h}}

+∥u^{h,2}∥_{Z^{h}} + ∥d_{t}u^{h,0}∥_{0,Ω_{p}} + ∥d_{t}u^{h,1}∥_{0,Ω_{p}} + ∑_{n=1}^{L−1} ∥F^{n}∥_{0,Ω_{p}}^{2}Δt) .

As a consequence of the estimate (132) we have uniqueness and existence for the procedure (113)-(115), provided the time step Δt satisfies the stability conditions (125) and (128).

7 The triangular case for the PSVTM-mode

Let T^{h} be a quasiregular partition of Ω into triangles Ω_{j} of diameter bounded by h. Set

R_{1} = {U : u = (a + b x_{2}, c + d x_{1}), a, b, cconstants},

and let us define the spaces to approximate the electric and magnetic fields V^{h} and V^{h} as

V^{h} = {ψ ∈ H(curl, Ω) : ψ|_{Ω_{j}} ∈ R_{1}(Ω_{j}) ∀j},
W^{h} = {φ ∈ L^{2}(Ω) : φ|_{Ω_{j}} ∈ P_{0}(Ω_{j}) ∀j}.

The local degrees of fredom for V^{h} and W^{h} can be taken to be the values of the tangential components at the mid-points of each edge of the triangle Ω_{j} and the values at the centroids of Ω_{j}, respectively.

Next, if

N C_{j}^{h} = P_{1}(Ω)_{j}),

the space N C^{h} to approximate each component of the solid displacement inΩ_{p} is defined as in the rectangular case as

N C^{h} = {v : v_{j} = v|_{Ω_{j}} ∈ N C_{j}^{h}, v_{j}(ξ_{jk}) = v_{k}(ξ_{jk}) ∀(j, k)},

with the local degrees of freedom being the values at the mid-point of the edges of each Ω_{j}.

Finally if

S_{1} = {U : u = (a + b x_{1}, c + d x_{2}), a, b, cconstants},

the space to approximate the fluid displacement vector is

M^{h} = {v^{f} ∈ H(div, Ω_{p}) : v^{f}|_{Ω_{j}} ∈ S_{1}(Ω_{j}) ∀j}.

The local degrees of fredom for M^{h} are the values of normal components at the mid-points of each edge of the triangle Ω_{j}.

The projection operators Π_{h}, P_{h}, Q_{h}, R_{h} and S_{h} are defined as in the rectangular case, with identical properties.

Next, setting

Y^{h} = V^{h} × W^{h} × (N C^{h})^{2} × M^{h}

the definition of the continuous and discrete-time finite element procedures remain unchanged, with identical apriori error estimates and stability results.

8 3D seismoelectric modeling. 3D rectangular elements.

Consider a 3D-rectangular domain Ω = Ω_{a}∪Ω_{p} where Ω_{a} and Ω_{p} are associated with the air and subsurface poroelastic (disjoint) parts of Ω, respectively. Consider the solution of

ε\frac{∂E}{∂t} + σE − ∇ × H + L_{0} \frac{η}{κ_{0}} \frac{∂u^{f}}{∂t} = 0, Ω,  (133)
∇ × E + \frac{∂H}{∂t} = 0, Ω,  (134)
ρ_{b} \frac{∂^{2}u^{s}}{∂t^{2}} + ρ_{f} \frac{∂^{2}u^{f}}{∂t^{2}} − ∇ · τ(u) = F^{(s)}, Ω_{p},  (135)
ρ_{f} \frac{∂^{2}u^{s}}{∂t^{2}} + m\frac{∂^{2}u^{f}}{∂t^{2}} + \frac{η}{κ_{0}} \frac{∂u^{f}}{∂t} u^{f} + ∇p_{f} = F^{(f)}, Ω_{p},  (136)

where τ(u), p_{f}(u) are given by (5) and (6), respectively. We will solve (133)- (136) with the boundary conditions

ε^{1/2}P_{χ}E + μ^{1/2}ν × H = 0, on Γ,  (137)
−G_{Γ_{p}}(u) = DS_{Γ_{p}}(\frac{∂u}{∂t}), on Γ_{p},  (138)

the free surface condition

−G_{Γ_{p}}(u) = 0, on Γ_{a,p},  (139)

and the initial conditions

E(t = 0) = E_{0}, H(t = 0) = H_{0}  (140)

u^{s}(t = 0) = u_{0}^{s}, u^{f}(t = 0) = u_{0}^{f},

\frac{∂u^{s}}{∂t} (t = 0) = u_{1}^{s}, \frac{∂u^{f}}{∂t} (, t = 0) = u_{1}^{f}.

Here

P_{χ}E = E − ν(ν · E) = −ν × (ν × E)

is the 3-D orthogonal projection of the trace of E into the tangential plane perpendicular to the normal vector ν and

G_{Γ_{s}}(u) = (τ(u)ν · ν, τ(u)ν · χ^{1}, τ(u)ν · χ^{2}, p_{f}(u))^{t},
S_{Γ_{s}}(u) = (u^{s} · ν, u^{s} · χ^{1}, u^{s} · χ^{1}, u^{f} · ν)^{t} ,

where , where ^{t} denotes the transpose, ν is the unit outer normal on Γ_{s} and χ^{1}, χ^{2} are two unit tangents on Γ_{s} such that {ν, χ^{1}, χ^{2}} is an orthonormal set on Γ_{p}. The matrices 4 × 4 matrices R, M in the definition of the matrix D in (138) are defined as the obvious extension to the 3D case of the defining equations in the PSVTM case (simply add a row and a column containing all zeroes and either b or G in the main diagonal as needed).

9 A weak formulation

First recall the integration by parts formula [42]

(∇ × U, V) − (U, ∇ × V) = ⟨ν × U, V ⟩ , ∀ U, V ∈ H(curl, Ω). (141)

Thus, multiply (133) by ψ ∈ H(curl, Ω), integrate in Ω and use (141). Also multiply (134) by φ ∈ [L^{2}(Ω)]^{3} and integrate. Finally, multiply (135) by v^{s} ∈ [H^{1}(Ω)]^{3} and (136) by v^{f} ∈ H(div, Ω), add the resulting equations, integrate in Ω_{p} and use (138). WE obtain the weak formulation: find (E, H, u^{s}, u^{f}) ∈ H(curl, Ω) × [L^{2}(Ω)]^{3} × [H^{1}(Ω)]^{3} × H(div, Ω) such that

(ε\frac{∂E}{∂t} , ψ) + (σE, ψ) − (H, ∇ × ψ) + (L_{0} \frac{η}{κ_{0}} \frac{∂u^{f}}{∂t} , ψ)_{Ω_{p}}  (142)

+ ⟨(\frac{ε}{μ})^{1/2} P_{χ}E, P_{χ}ψ⟩ = 0, ψ ∈ H(curl, Ω),

(∇ × E, φ) + (μ\frac{∂H}{∂t} , φ) = 0, φ ∈ [L^{2}(Ω)]^{3}.  (143)
(P \frac{∂^{2}u}{∂t^{2}}, v)_{Ω_{p}} + (\frac{η}{κ_{0}} \frac{∂u^{f}}{∂t} , v^{f})_{Ω_{p}} + A(u, v) + ⟨DS_{Γ_{p}}(\frac{∂u}{∂t}), S_{Γ_{p}}(v)⟩_{Γ_{p}} (144)
= (F, v)_{Ω_{p}}, v = (v^{s}, v^{f}) ∈ [H^{1}(Ω_{p})]^{3} × H(div, Ω_{p}).

In (144) the matrix P is defined as in (24) and A(u, v) is defined as in (26) provided the symmetric 7 × 7 matrix M = (m_{ij}) is defined with entries m_{11} = m_{22} = m_{33} = λ_{c} + 2G, m_{12} = m_{13} = λ_{c}, m_{14} = αK_{av}, m_{15} = m_{16} = m_{17} = 0, m_{23} = λ_{c}, m_{24} = αK_{av}, m_{25} = m_{26} = m_{27} = 0, m34 = αK_{av}, m_{35} = m_{36} = m_{37} = 0, m_{44} = K_{av}, m_{45} = m_{46} = m_{47} = 0, m_{55} = m_{66} = m_{77} = G, m_{56} = m_{57} = m_{67} = 0. Also, define

ε̃(u) = (ε_{11}(u^{(s)}), ε_{22}(u^{(s)}), ε_{33}(u^{(s)}), ∇ · u^{(f)}, ε_{12}(u^{(s)}, ε_{13}(u^{(s)}, ε_{23}(u^{(s)})^{t}.

Now equations (142)-(144) are formally identical to equations (22)-(25) for the 2D case and the argument for uniqueness and existence of the continuous problem can be repeated to demonstrate that the statement in Theorem 1 remains valid in the 3D case.

10 A continuous-time finite element procedure

Let T^{h} be a quasiregular partition of Ω into parallelepipeds Ω_{j} of diameter bounded by h. To approximate the electromagnetic fields E, H we will employ the 3D mixed finite element space V^{h} × W^{h}, defined as follows [18,21]:

V^{h} = {ψ ∈ H(curl, Ω) : ψ|Ω_{j} ∈ V_{j}^{h} ≡ P_{0,1,1} × P_{1,0,1} × P_{1,1,0}},
W^{h} = {φ ∈ L^{2}(Ω) : φ|Ω_{j} ∈ W_{j}^{h} ≡ P_{1,0,0} × P_{0,1,0} × P_{0,0,1}}.

Note that functions in V^{h} have continuous projections of their traces across the interelemnt boundaries Γ_{jk}. Also,

∇ × V^{h} ⊂ W^{h}.

For any ψ ∈ V^{h} and any φ ∈ W^{h} the local degrees of freedom in Ω_{j} are [21]

∑_{j}^{ψ} = {(ψ · χ_{p}^{j})(m_{p}), where m_{p} is the mid-point of the p-edge  (145)
e_{p}^{j} of Ω_{j} of unit tangent χ_{p}^{j}, 1 ≤ p ≤ 12} ,
∑_{j}^{φ} = {(φ · ν_{p}^{j})(G_{p}), where G_{p} is the centroid of the p-face  (146)
f_{p}^{j} of Ω_{j} of unit normal ν_{p}^{j}, 1 ≤ p ≤ 6} .

This definition of the spaces V^{h}, W^{h} and their degrees of freedom is in pages 114-115 of [21], formulas (60), 63), (66) and (67). Monk calls U_{h} to V^{h} and V̂_{h} to W^{h}.

Let us define the projection

Π_{h} : H(curl, Ω) → V^{h} : ∫_{e_{j}^{p}}(ψ − Π_{h}ψ) · χ_{p}^{j}ds = 0, 1 ≤ p ≤ 12,  (147)

(148)

for each p-edge e_{p}^{j} with tangent χ_{p}^{j} of Ω_{j}.

Also, let the L^{2}-projection be defined as in the 2D-case in (41). Then the approximating properties (45), (46) and (47) remain valid.

The nonconforming finite element space to approximate the solid displacement vector is defined as follows [23]. On the reference element R̂ = [−1, 1]^{3} set

Q(R̂) = Span\{1, x̂_{1}, x̂_{2}, x̂_{3}α̃(x̂_{1}) − α̃(x̂_{2}), α̃(x̂_{1}) − α̃(x̂_{3})\}, α̃(ˆ1) = x̂_{1}^{2} − \frac{5}{3}x̂_{1}^{4}.

Then let N C_{j}^{h} = Q(Ω_{j}) and

N C^{h} = {v : v_{j} = v|_{Ω_{j}} ∈ N C_{j}^{h}, v_{j}(ξ_{jk}) = v_{k}(ξ_{jk}) ∀(j, k)}.

The six local degres of freedom are the values at the centers ξ_{jk} of the faces of Ω_{j}.

To approximate the fluid displacement, we employ the vector part of the Raviart-Thomas-Nedelec space of zero order [17,18], defined as follows.

M^{h} = {w ∈ H(div, Ω_{p}) : w|Ω_{j} ∈ M_{j}^{h} ≡ P_{1,0,0} × P_{0,1,0} × P_{0,0,1}}.

The projections

R_{h} : [H^{2}(Ω_{p})]^{3} → [N C^{h}]^{3},
Q_{h} : [H^{1}(Ω_{p})]^{3} → M^{h},
S_{h} : [H^{2}(Ω_{p})]^{3} × H^{1}(div; Ω_{p}) → Λ̃^{h},

are defined identically than in (42), (43) and (44), where now Λ̃^{h} is

Λ̃^{h} = {λ̃^{h} : λ̃_{jk}^{h} = tr _{Γ_{jk,p}}(λ̃^{h}|_{Ω_{j}}) ∈ [P_{0}(Γ_{jk,p})]^{3} ≡ Λ̃_{jk}^{h}, λ̃_{jk}^{h} + λ̃_{kj}^{h} = 0} ,

the approximatting properties (48), (49) and (50) still hold for these spaces. Set

Θ_{h} ((E, H, u^{s}, u^{f}), (ψ, φ, v^{s}, v^{f})) = (ε\frac{∂E}{∂t} , ψ) + (σE^{,}ψ) − (H, ∇ × ψ)  (149)

+ (L_{0} \frac{η}{κ_{0}} \frac{∂u^{f}}{∂t} , ψ)_{Ω_{p}} + (∇ × E, φ) + ⟨(\frac{ε}{μ})^{1/2} E · χ, ψ · χ⟩

+(μ\frac{∂H}{∂t} , φ) + (P \frac{∂^{2}u}{∂t^{2}}, v)_{Ω_{p}} + (\frac{η}{κ_{0}} \frac{∂u^{f}}{∂t} , v^{f} )_{Ω_{p}} + A_{h}(u, v) + ⟨DS_{Γ}(\frac{∂u}{∂t}), S_{Γ}(v)⟩_{Γ_{p}} ,

where A_{h}(u, v) is defined as in (53) using the 3D definition given here for the matrix M.

Setting

Y^{h} = V^{h} × W^{h} × (N C^{h})^{3} × M^{h},

the definition of the continuous-time Galerkin procedure (55)-(57) remains unchanged, as well as the existence and uniquess results and apriori error estimates in Theorem 2 and Theorem 3.

11 The discrete-time finite element procedures for 3D rectangular elements

The first fully implicit discrete-time procedure is the analogue of (101)-(103) and is formulated as follows. Given (E^{h,1}, H^{h,1}, u^{s,h,0}, u^{s,h,1}, u^{f,h,0}, u^{f,h,1}) ∈

V^{h}×W^{h}×(N C^{h})^{6}×(M^{h})2, for n ≥ 1 compute (E^{h,n+1}, H^{h,n+1}, u^{s,h,n+1}, u^{s,h,n+1}) ∈ Y^{h} such that

(ε d_{t}E^{h,n}, ψ) + (σÊ^{h,n+1/2}, ψ) − (Ĥ^{h,n+1/2}, ∇ × ψ)  (150)

+ (L_{0} \frac{η}{κ_{0}} ∂u^{f,h,n}, ψ)_{Ω_{p}} + ⟨(\frac{ε}{μ})^{1/2} P_{χ}Ê^{h,n+1/2}, P_{χ}ψ⟩ = 0, ψ ∈ V^{h},

(∇ × Ê^{h,n+1/2}, φ) + (μd_{t}H^{h,n}, φ) = 0, φ ∈ W^{h},  (151)
(P∂^{2}u^{h,n}, v)_{Ω_{p}} + (\frac{η}{κ_{0}} ∂u^{f,h,n}, v^{f} )_{Ω_{p}} + A_{h}(\frac{u^{h,n+1} + u^{h,n−1}}{2} , v)  (152)
+ ⟨DS_{Γ}(∂u^{h,n}), S_{Γ}(v)⟩_{Γ_{p}} = (F^{n}, v), v = (v^{s}, v^{f}) ∈ (N C^{h})^{3} × M^{h}.

The other two discrete-time procedures can be formulated with the obvious changes for the 3D case, as it is done in (150)-(152). The analysis on existence, uniqueness and stability follows for these 3D procedures follow with identical arguments than for the 2D case and is omitted for brevity.

12 The case of tetrahedral elements

Let T^{h}(Ω) be a nonoverlapping quasiregular partition of Ω = Ω_{p} ∪ Ω_{a} into tetrahedral elements Ω_{j} of diameter bounded by h and for x = (x_{1}, x_{2}, x_{3}) let

R_{1} = {U : u = (α + β x, α, β ∈ (P_{0})^{3}},

Then following [19], [20] the space V^{h} to approximate the electric field V^{h} is

V^{h} = {ψ ∈ H(curl, Ω) : ψ|_{Ω_{j}} ∈ R_{1}(Ω_{j}) ∀j}

The space to approximate the magnetic field is

W^{h} = {φ ∈ (L^{2}(Ω))^{3} : φ|_{Ω_{j}} ∈ (P_{0}(Ω_{j}))^{3} ∀j}

Nedelec [19] showed that the following degrees of freedom

M_{j}^{ψ} = {∫_{e_{j}} ψ · χ^{j}ds, for the six edges e_{j} of Ω_{j}} ,

where χ_{j} is a unit vector parallel to e_{j} are R_{1}-solvent and curl-conforming, that curl V^{h} ⊂ W^{h} and that (46) holds.

The space V^{h} and W^{h} and their DOF defined here are defined as here in [20] in formulas (3.2), (3.3) of page 1618 and (3.8) of page 1620. They are also defined in [21] in pages 108-109. See what he says about the choice of his V_{h} that is our W^{h}, he says that is different from the Nedelec space, and that his choice, that is our choice, is better because the matrix is block diagonal.

Agregar la definicion de los DOf para el espacio W^{h}, debe ser una pelotudez pero no esta escrito en nungun paper de los que refiero....

Let us define the projection

Π_{h} : H(curl, Ω) → V^{h} : ∫_{e_{j}}(ψ − Π_{h}ψ) · χ^{j}ds = 0, 1 ≤ j ≤ 6,  (153)

(154) where χ^{j} is a unit vector parallel to e_{j}.

Also, let the L^{2}-projection be defined as in the 2D-case in (41). Then the approximating properties (45), (46) and (47) remain valid.

The nonconforming finite element space to approximate the solid displacement vector is defined as in [23]. Let N C_{j}^{h} = P_{1}(Ω_{j}) and

N C^{h} = {v : v_{j} = v|_{Ω_{j}} ∈ N C_{j}^{h}, v_{j}(ξ_{jk}) = v_{k}(ξ_{jk}) ∀(j, k)}.

The four local degres of freedom are the values at the centers ξ_{jk} of the faces of Ω_{j}. Next, let

S_{1} = {U : U = α + λx, α ∈ (p_{0})^{3}, λ ∈ P_{0}}  (155)

and let

M^{h} = {v^{f} ∈ H(div, Ω_{p}) : v_{j} = v|_{Ω_{j}} ∈ S_{1} ∀j}.  (156)

The projections

R_{h} : [H^{2}(Ω_{p})]^{3} → [N C^{h}]^{3},
Q_{h} : [H^{1}(Ω_{p})]^{3} ∪ (H_{h}^{1}(Ω_{p}))^{3} → M^{h},
S_{h} : [H^{2}(Ω_{p})]^{3} × H^{1}(div; Ω_{p}) → Λ̃^{h},

are defined identically than in (42), (43) and (44). Nedelec [19] showed that the degrees of freedom

N_{j}(v^{f}) = {⟨v^{f} · ν, 1⟩_{B} , B any of the four faces of Ω_{j}}  (157)

are S_{1}-solvent and conforming in H(div; Ω_{p}), so that (43) uniquely defines Q_{h}.

Setting

Y^{h} = V^{h} × W^{h} × (N C^{h})^{3} × M^{h},

the definition of the continuous-time Galerkin procedure (55)-(57) remains unchanged, as well as the existence and uniquess results and apriori error estimates in Theorem 2 and Theorem 3.

13 3D seismoelectric modeling

Creo que aqui debemos definir el modelo 3D para SE con una fuente J_{s}^{e} ∈ L^{2} in formular la forma debil y demostrar unicidad usando algo analogo al caso SE. Tambien definir los metodos de elementos finitos continuos y discretos y decir que el analisis es analogo.

Despues en otro paper nos ocupamos de como se ponen los electrodos, etc.

14 Numerical experiments for 2D seismoelectric modeling

Aca se tiene que hamacar fabio y Patricia para sacar rapido este paper.

Miren el apper de Cohen y Monk en CMAME, vol 169, 1999, p197-217, alli en la seccion 3.3 usa Leapfrog para Maxwell. tambien usan Leapfrog y dan una cond. de estabilidad en el paper de cohen y monk que les mando. Parece que hay que usar Lepafrog.

References

[1] S. R. Pride and M. W. Haartsen, “Electroseismic wave properties,” J. Acoust. Soc. Amer., 100 (3), 1301 (1996).

[2] O. V. Mikhailov, M. W. Haartsen and M. N. Toksoz, “Electroseismic investigation of the sallow subsurface: Field measurements and numerical modeling,” Geophysics, 62, 97 (1997).

[3] O. V. Mikhailov, J. Queen and M. N. Toksoz, “Using borehole electroseismic measurements to detect and characterize fractured (permeable) zones,” Geophysics, 65, 1098 (2000).

[4] A. H. Thompson, “Electromagnetic-to-seismic conversion: Successful developments suggest viable applications in exploration and production,” 75th SEG Annual Meeting, Expanded Abstracts, Houston, (2005).

[5] S. Hornbostel and A. H. Thompson, “Waveform design for electroseismic exploration,” Geophysics, 72 (2), Q1 (2007)

[6] A. H. Thompson and G. Gist, “Geophysical applications of electrokinetic conversion,” The Leading Edge, 12, 1169 (1993).

[7] S. R. Pride, “Governing equations for the coupled electromagnetics and acoustics of porous media,” Physics Review B, 50, 15678 (1994), [8] G. I. Block and J. G. Harris, “Conductivity dependence of seismoelectric wave

phenomena in fluid-saturated sediments,” J. Geophys. Res., 111, B01304 (2006), DOI:10.1029/2005JB003798.

[9] S. H. Haines and S. R. Pride, “Seismoelectric numerical modeling on a grid,” Geophysics, 71 (6), N57 (2006)

[10] M. A. Biot, “Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low frequency range,” J. Acoust. Soc. Amer., 28, 168 (1956).

[11] M. A. Biot, “Theory of propagation of elastic waves in a fluid-saturated porous solid. II. High frequency range,” J. Acoust. Soc. Amer., 28, 179 (1956).

[12] M. H. Haartsen and S. R. Pride, “Electroseismic waves from point sources in layered media,” J. Geophys. Res., 102, (NO. B11), 24745 (1997).

[13] Q. Han and Z. Wang, “Time-domain simulation of SH-wave-induced electromagnetic field in heterogeneous porous media: A fast finite element algorithm,” Geophysics, 66 (2), 448 (2001)

[14] C. C. Pain, J. H. Saunders, M. H. Wortington, J. M. Singer, W. Stuart-Bruges, G. Mason and A. Goddard, “A mixed finite element method for solving the poroelastic Biot equations with electrokinetic coupling,” Geophys. J. Int., 160, 592 (2005).

[15] B. S. White, “Asymptotic theory of electroseismic prospecting,” SIAM J. Appl. Math., 65 (4), 1443 (2005).

[16] B. S. White and M. Zhou, “Electroseismic prospecting in layered media,” SIAM J. Appl. Math., 67 (1), 69 (2006).

[17] P. A. Raviart and J. M. Thomas, “A mixed finite element method for second order elliptic problems,” in Mathematical Aspects of the Finite Element Method, Lecture Notes in Mathematics, volume 606, Springer-Verlag, Berlin, New York, I. Galligani and E. Magenes, eds., (1977), p.292.

[18] J. C. Nedelec, “Mixed finite elements in R^{3},” Numer. Math., 35, 315 (1980).

[19] J. C. Nedelec, “A new family of mixed finite elements in R^{3},” Numer. Math.,50, 57 (1986).

[20] P. B. Monk, “A mixed method for approximating Maxwell’s equations” SIAM J. Num. Anal., 28 (6), 1610 (1991).

[21] P. B. Monk, “An analysis of Nedelec method for the spatial discretization of Maxwell equations” J. Comput. Appl. Math, 47, 101 (1993).

[22] P. B. Monk, and A. K. Parrot, “A dispersion analysis of finite element methods for Maxwell’s equations,” SIAM J. Sci. Stat. Comput., 15 (4), 916 (1994).

[23] J. Douglas, Jr., J. E. Santos, D. Sheen and X. Ye, “ Nonconforming Galerkin methods based on quadrilateral elements for second order elliptic problems,” RAIRO Mathematical Modelling and Numerical Analysis (M2AN) 33, 747 (1999).

[24] J. E. Santos, “Elastic wave propagation in fluid-saturated porous media. Part I. The existence and uniqueness results,” Mathematical Modelling and Numerical Analysis (M2AN) , 20 (1), 113 (1986).

[25] J. L. Lions, “ Quelques methodes dde resolution des problems aux limites nonlineaires”, Dunod, Gauthier-Villars, Paris, (1969).

[26] J. E. Santos and D. Sheen, “Finite element methods for the simulation of waves in composite saturated poroviscoelastic materials,” SIAM J. Numer. Anal., 45 (1), 389 (2007).

[27] F. I. Zyserman and J. E. Santos, “Analysis of the numerical dispersion of waves in saturated poroelastic media,” Computer Methods in Applied Mechanics and Engineering, 196, 4644 (2007).

[28] J. Douglas, Jr., J. E. Santos, and D. Sheen, “A nonconforming mixed finite element method for Maxwell’s equations,” Mathematical Models and Methods in Applied Sciences, 10, 593 (2000).

[29] J. E. Santos, “Global and domain-decomposed mixed methods for the solution of Maxwell’s equations with application to magnetotellurics,” Numer. Methods Partial Differ. Eqs., 14, 1 (1998).

[30] J. E. Santos and D. Sheen, “On the existence and uniqueness of Maxwell’s equations in bounded domains with application to magnetotellurics,” Mathematical Models and Methods in Applied Sciences, 10, 615 (2000).

[31] F. I. Zyserman, L. Guarracino and J. E. Santos, “A hybridized mixed finite element domain decomposed method for two-dimensional magnetotelluric modelling,” Earth, Planets and Space, 51, 297 (1999).

[32] F. I. Zyserman and J. E. Santos, “Parallel finite element algorithm for three– dimensional magnetotelluric modelling,” J. Appl. Geophysics, 44, 337 (2000).

[33] D. L. Johnston, J. Koplik and R. Dashen, “Theory of dynamic permeability and tortuosity in fluid–saturated porous media,” J. Fluid Mechanics, 176, 379 (1987).

[34] F. Gassmann, “Uber die elastizität poröser medien,” Vierteljahrsschrift der Naturforschenden Gessellschaft in Zurich, 96, 1 (1951).

[35] J. E. Santos, J. Douglas Jr, M. E. Morley and O. M. Lovera, “Reflection and transmission coefficients in fluid-saturated porous media,” J. Acoust. Soc. Amer., 91, 1911 (1992).

[36] M. A. Biot, “Theory of deformation of a porous viscoelastic anisotropic solid,” J. Appl. Physics 27, 459 (1956).

[37] M. A. Biot, “Mechanics of deformation and acoustic propagation in porous media,” J. Appl. Physics 33 (4), 1482 (1962).

[38] H. P. Liu and D. L. Anderson and H. Kanamori “Velocity dispersion due to anelasticity; implications for seismology and mantle composition,” Geophys. J. R. Astr. Soc. l47, 41 (1976).

[39] S. Ward and G. W. Hohmann, Electromagnetic theory for geophysical applications, in Electromagnetic Methods in Applied Geophysics, Vol. 1, Theory, M. N. Nabiguian, ed., SEG Investigations in Geophysics # 3, (1987).

[40] J. E. Santos, J. Douglas, Jr., M. E. Morley, and O. M. Lovera, “Finite element methods for a model for full waveform acoustic logging”, IMA Journal of Numerical Analysis, 8, 415 (1998).

[41] V. Girault and P. Raviart, Finite Element Methods for Navier-Stokes Equations, Springer–Verlag, Berlin, 1986.

[42] D. Sheen, “A generalized Green’s theorem,” Appl. Math. Lett., 5, 95 (1992).

[43] D. Sheen, “Approximation of electromagnetic fileds: Part I. Continuous problems,” SIAM J. Appl. Math., 6, 1716-1736 (1997).

[44] A. Buffa and P. Ciarlet, “On traces for functional spaces related to maxwell’s equations. II. Hodge decompositions on the boundary of Lipschitiz polyhedra and applications,” Math. Methods Appl. Sci., 24 (1), 31 (2001).

[45] A. Buffa and M. Costabel and D. Sheen, “On traces for H(curl, Ω) in lipschitz domains,” J. math. Anal. Appl., 276 (2), 845 (2002).

[46] C. L. Ravazzoli and J. E. Santos, “A theory for wave propagation in porous rocks saturated by two-phase fluids under variable pressure conditions,” Bollettino di Geofisica teorica ed applicata, 46 (4), 261 (2005).

[47] , P. G. Ciarlet, The Finite Element Method for Elliptic Equations , North– Holland, Amsterdam , 1978.

[48] R. Durán, “Finite elements in H(curl) with applications to elasticity,” Numer. Methods Partial Differ. Eqs., 4, 167 (1990).

[49] J. Douglas Jr., and J. E. Roberts, “Global estimates for mixed methods for second order elliptic equations,” Math. Comput., 44, 39 (1985).

[50] Thomas, J. M., “Sur l’analyse numérique des methodes d’éléments finis hybrides et mixtes,” Ph.D. thesis, Université Pierre-et-Marie Curie, Paris, (1977).

[51] G. Duvaut and J. L. Lions, Les Inéqualities en Méchanique et en Physique, Dunod, Paris, 1972.

[52] J. A. Nitsche, “On Korn’s second inequality,” R. A. I. R. O Anal. Numer., 15, 237 (1981).