Modeling and simulation of nonlinear electro-thermo-mechanical continua with application to shape memory polymeric medical devices
MModeling and simulation of nonlinearelectro-thermo-mechanical continua with application toshape memory polymeric medical devices
Innocent Niyonzima a , Yang Jiao a , Jacob Fish a a Columbia University, Department of Civil Engineering and Engineering Mechanics, NewYork, 10027 NY, USA.
Abstract
Shape memory materials have gained considerable attention thanks to theirability to change physical properties when subjected to external stimuli such astemperature, pH, humidity, electromagnetic fields, etc. These materials are in-creasingly used for a large number of biomedical applications. For applicationsinside the human body, contactless control can be achieved by the addition ofelectric and/or magnetic particles that can react to electromagnetic fields, thusleading to a composite biomaterial. The difficulty of developing accurate nu-merical models for smart materials results from their multiscale nature and fromthe multiphysics coupling of involved phenomena. This coupling involves elec-tromagnetic, thermal and mechanical problems. This paper contributes to themultiphysics modeling of a shape memory polymer material used as a medicalstent. The stent is excited by electromagnetic fields produced by a coil whichcan be wrapped around a failing organ. In this paper we develop large defor-mation formulations for the coupled electro-thermo-mechanical problem usingthe electric potential to solve the electric problem. The formulations are thendiscretized and solved using the finite element method. Results are validatedby comparison with results in the literature.
Keywords:
Multiphysics modeling, electro-thermo-mechanical coupling, shapememory polymers stents, large deformations.
Preprint submitted to Computer Methods in Applied Mechanics and Engineering April 29, 2019 a r X i v : . [ m a t h . NA ] A p r . Introduction The increase of life expectancy creates a need to maintain the functions of aging organs to allow greater independence for the elderly. Biomaterial implants have the potential to fulfill some of these functions. The total number of im- plants in the world exceeds four hundred million per year and grows every year [1]. Biomaterials are also increasingly used for a large number of biomedical ap- plications such as the prevention and cure of coronary heart disease and stroke, as well as ophthalmological applications, biosensors and drug delivery systems [2, 3]. They have the potential to contribute to the reduction of the cost of health and the improvment of the life conditions. Among biomaterials, shape memory materials have gained considerable at- tention in the biomedical community thanks to their ability to change physical properties (morphing, structural rigidity, refractive index, etc.) when subjected to external stimuli such as temperature, pH, humidity, electromagnetic fields, etc. This special behavior results from the the shape memory effect observed in shape memory materials [4]. They are used in minimally invasive surgery as embolic devices to treat aneurysm [5, 6] and as vascular stents [7, 8, 9]. Fig- ure 1 illustrate the deployment of a stent in a blood vessel. They can also be used as portable sensors to monitor heart and respiratory rates and in con- trolled drug delivery systems thus allowing to reduce the side effects of drugs [11, 12]. For these different uses, biomaterials must possess a number of prop- erties. They must be biocompatible to avoid toxicity in contact with biological tissues. Biodegradability is a desirable property for temporary implants, and for minimally invasive in vivo surgery applications, devices must be controllable without contact and self-expanding. All of these properties make polymers the best candidates for a wide range of biomedical applications. Contactless control can be achieved by the addition of electric/magnetic (nano)particles inclusions to produce a smart composite that can react to electromagnetic fields. Accurate numerical models for smart composites must account for the multi- physics coupling which involve different domains of physics (electromagnetism, igure 1: Diagram illustrating the deployment of a stent in a blood vessel [10]. thermal, mechanics) and the multiscale nature of the materials. The present manuscript is concerned with multiphysics modeling of shape memory polymer materials used for biomedical devices for homogeneous materials. The difficulty of developing multiphysics models arises from a number of factors. (a) The geometric non-linearity resulting from large mechanical deformations leads to the modification of the equations that govern the electromagnetic and thermal problems to account for motion. Additional complexity for the electromagnetic problem results from the presence of electromagnetic fields in the air and vac- uum. (b) The material non-linearities resulting from the presence of materials with nonlinear thermal constitutive laws, plastic and viscous laws for mechanics and nonlinear anhysteretic/hysteretic constitutive laws for electromagnetism. Models and numerical simulations have already been developed for mul- tiphysics problems in piezoelectric, magnetostrictive, and piezomagnetic ma- terials in the case of small deformations [13, 14, 15, 16, 17, 18, 19]. The- oretical models have also been developed for thermomechanical and electro- magneto-thermomechanical problems in the case of large mechanical deforma- grangian approach for the mechanical problem and the Eulerian or arbitrary Lagrangian Eulerian approach for the electromagnetic problem [24, 25]. In this paper we consider low frequency electromagnetic problems and solve for a scalar potential formulation only defined in the mechanical domain. Thus, a Lagrangian mesh can also be used for the electromagnetic problem. The development of multiphysics models for smart composites controlled by electromagnetic fields is still in its infancy. Electro-magneto-mechanical mod- els using electrostatic and magnostatic formulations have been developed for magneto-sensitive composites and magneto-electro-elastic composites (e.g., elec- tro active polymers ) in large deformations [26, 27, 28]. In [29], the discontinuous Galerkin method was used to solve the electro-thermo-mechanical problem in a SMP by solving an electrokinetic problem excited by surface currents. For a more effective contactless control, the multiphysics problem should include eddy currents and hysteretic losses as a means of controlling the temperature. In this paper, we develop a simple multiphysics model for an electromag- netically controlled vascular stent excited by a coil. The paper extends the thermomechanical model developed in [30] by proposing a contactless electro- magnetic control of the temperature, especially during the recovery step that takes place inside the human body. For the sake of clarity and in order to have a self-sufficient paper which is easily accessible by the mechanics and electromag- netic communities, we derive the general fully coupled problem from Maxwell equations and conservation laws using the Lagrangian and Eulerian formalisms. Then, a simplified, quasistatic electro-thermo-mechanical problem is derived and discretized using the finite element (FE) method. The paper is organized as follows: in Section 2 we recall Maxwell’s equations and conservation laws using the Lagrangian and the Eulerian formalisms. In Sec- tion 3, we derive the simplified coupled problem and its strong and weak forms using potential formulations. The weak formulations are then semi-discretized in space using the FE method and in time using the backward Euler time step- ping method. The resulting system of nonlinear algebraic equations is linearized cal examples. At first, we validate the thermo-mechanical formulation for SMP materials with simple geometries along the lines of [30]. We then study the behavior of the electromagnetically responsive SMP stent excited by a coil. In Section 5 we close the paper with conclusions and perspectives.
2. Governing equations of the general multiphysics problem In this section, the general electro-magneto-thermo-mechanical coupled prob- lem is derived from Maxwell’s equations and conservation laws. Throughout the paper, we use the indices E and L to denote the Eulerian and Lagrangian quan- tities. Thus, f E and f L denote the forces in the Eulerian and Lagrangian framework, respectively. The open domains Ω Mec0 , Ω The0 and Ω Ele0 denote the undeformed computational domains for the mechanical, thermal and electro- magnetic problems, respectively. Likewise, Ω Mect , Ω Thet and Ω Elet denote the deformed computational domains for the mechanical, thermal and electromag- netic problems at a time t ∈ I t :=] t , t end [ . The domains of the mechanical and thermal problems are generally subdomains of the electromagnetic domain, i.e., Ω Mec i ⊆ Ω Ele i and Ω The i ⊆ Ω Ele i with i = { , t } as the electromagnetic fields can be defined in the entire domain including the surrounding air. The domains Ω c,i ⊆ Ω Ele i , Ω Cc,i ⊆ Ω Ele i and Ω s,i (cid:40) Ω Cc,i ⊂ Ω Ele i are the conductors, non- conductors and inductors where the electric currents source is imposed. The domains Γ Ele i , Γ Ele i and Γ Ele i denote the boundaries of the electromagnetic, ther- mal and mechanical domains, respectively. The differential operators Grad , Curl and Div denote the gradient, rotational and divergence operators defined on the undeformed configurations while grad , curl and div denote the same operators defined on the deformed configurations. The motion is described by the mappings ϕ t and ϕ assumed to be smoothenough (we do not consider fracture). The mapping ϕ t is also assumed to be5ijective and defined by: ϕ t : Ω Mec0 → E , X (cid:55)→ x = ϕ t ( X ) = ϕ ( X , t ) = X + u ( X , t ) (2.1)where X is the position of a particle point P in the undeformed configuration, x is the position of P in the deformed configuration, u is the vector of dis-placements and E is the three dimensional Euclidean space [31]. The positionsin the undeformed and deformed configurations are related by X = ϕ − t ( x ) which is valid thanks to the bijection of ϕ t . For any time t ∈ I t , the deformedconfigurations are also defined as: Ω Mect := ϕ t (cid:0) Ω Mec0 (cid:1) , Ω Thet := ϕ t (cid:0) Ω The0 (cid:1) , Ω Elet := ϕ t (cid:0) Ω Ele0 (cid:1) . (2.2)The deformation gradient tensor and its determinant are given by: F := ∂ x ∂ X = Grad ϕ = + Grad u , J = det F (2.3)where is the identity matrix. The velocity v and acceleration a are given by: v ( X , t ) = ∂ ϕ ∂t ( X , t ) , a ( X , t ) = ∂ v ∂t ( X , t ) = ∂ ϕ ∂t ( X , t ) . (2.4)Assuming the existence of a mapping Θ : Θ : E × I t → Ω Mec0 (cid:40) E , ( x , t ) (cid:55)→ X = Θ ( x , t ) = Θ ( ϕ ( X , t ) , t ) , (2.5)it is possible to derive the following relationship: D X Dt = ∂ Θ ∂ ϕ ∂ ϕ ∂t + ∂ Θ ∂t = F − v + V = , (2.6)which relates the matter flow field V to the velocity v as: V = − F − v . (2.7)The matter flow field is important for the definition of the electromagnetic prob- lem on the undeformed configuration. The independence of the initial position X on the time was used to derive (2.6). .2. Maxwell’s equations We recall non-relativistic Maxwell’s equations on moving domains under large deformations, using Eulerian and Lagrangian formalisms [32, 33, 21, 20].
In Eulerian setting, the electromagnetic fields are governed by the followingMaxwell’s equations [34, 20]: curl h = j + ∂ t d , curl e = − ∂ t b , div b = 0 , div d = ρ E , (2.8 a-d)and constitutive laws: h = µ − b − ( m − v × p ) (cid:124) (cid:123)(cid:122) (cid:125) m eff = ν E ( b ) b + v × p = H ( e , b , v ) , d = (cid:15) e + p = (cid:15) (cid:15) rE ( e ) e = (cid:15) E ( e ) e = D ( e ) , InnoInnoInnoInno j = σ E ( e + v × b (cid:124) (cid:123)(cid:122) (cid:125) e eff ) + j s + ρ E v = j eff + ρ E v = J ( e , b , v ) . (2.9 a-c)In (2.8 a-d)–(2.9 a-c), h is the magnetic field (A/m), b the magnetic flux density(T), j the electric current density (A/m ), d the electric flux density (C/m ), e the electric field (V/m) and ρ E the electric charge density (C/m ). The fields e eff and j eff are the effective electric field and effective electric current density,whereas j s (A/m ) is the electric current source defined in the inductors Ω s,t and σ E is the electric conductivity tensor ( Ω /m). The magnetization m (A/m)and polarization p (C/m ) are defined by: m = µ − χ bE ( b ) b , p = (cid:15) χ eE ( e ) e = χ eE + χ eE d , (2.10 a-b)where χ bE and χ eE are the magnetic and electric susceptibility tensors, µ = π − is the magnetic permeability of the free space (H/m) and (cid:15) (cid:39) − / π is the electric permittivity of the free space (C /Nm ). Another definition of magnetic susceptibility χ mE with m = χ mE h is often used in the constitu- tive law dual to (2.10 a). Additionally, ν E = µ − ( − χ bE ) is the magnetic reluctivity tensor, µ E = ν − E is the magnetic permeability tensor, (cid:15) E is the (cid:15) rE ( e ) is the relative electric permittivity tensor with (cid:15) rE ( e ) = + χ eE ( e ) and v is the velocity (m/s). The dependency of the mapping H on (2.9 a) on the electric field e results from (2.10 b). The motion is accounted for by the velocity terms in the constitutive laws (2.9 a-c). In Lagrangian setting, the electromagnetic fields are governed by the follow-ing Maxwell’s equations (see [32, 33, 21, 20, 23] and [35, Appendix F]):
Curl H eff = J + ∂ t D , Curl E eff = − ∂ t B , Div B = 0 , Div D = ρ L , (2.11 a-d)and constitutive laws: H eff = J − F T ( ν E ◦ ϕ − t ) F (cid:124) (cid:123)(cid:122) (cid:125) ν L B + ( (cid:15) rE ◦ ϕ − t ) − ( V × D ) (2.12) D = J F − ( (cid:15) E ◦ ϕ − t ) F − T (cid:124) (cid:123)(cid:122) (cid:125) (cid:15) L ( E eff + V × B (cid:124) (cid:123)(cid:122) (cid:125) E ) , (2.13) J = J F − ( σ E ◦ ϕ − t ) F − T (cid:124) (cid:123)(cid:122) (cid:125) σ L E eff + J s (2.14)defined on the undeformed configuration Ω Ele0 . In (2.11 a-d) and (2.12)–(2.14), H is the magnetic field (A/m), B the magnetic flux density (T), J the electric current density (A/m ), D the electric flux density (C/m ), E eff the effective electric field, E the electric field (V/m), J s (A/m ) the current source defined in the inductors Ω s, and V is the matter flow field (m/s) defined in (2.6). The tensor σ L is the electric conductivity tensor ( Ω /m) and ρ L is the electric charge density (C/m ). The notation ( f ◦ ϕ t − )( x ) := f ( ϕ t − ( x )) is used. The one differential forms are transformed as: H eff = F T ( h − v × d ) , E eff = F T ( e + v × b ) , E = E eff + V × B = F T e . (2.15 a-c)and the two differential forms are transformed as: B = J F − b , D = J F − d , J = J F − j , J s = J F − j s . (2.16)8n (2.15 c), V × B = − F T ( v × b ) results from the identity: F T v × F T b = J F − ( v × b ) , (2.17)which is valid for any matrix F ∈ GL ( R ) [35, Formula B.11]. Additionally, the magnetization M is related to the magnetic induction by: M = F T m = F T ( µ − χ bE b ) = F T ( µ − χ bE J − F B )= µ − J − F T χ bE F B = µ − χ bL B = ( µ − − ν L ) B , (2.18)while the effective magnetization M eff is given by: M eff = F T ( m − v × p ) = M + V × P , = ( µ − − ν L ) B + V × P (2.19)where the polarization is related to the electric flux density by [35]: P = (cid:15) rL − (cid:15) rL D = χ eE + χ eE D . (2.20)Combining all these results, the following transformations for second order ten-sors used in constitutive laws can be derived: ν L = J − F T ν E F , µ L = J F − µ E F − T , (cid:15) L = J F − (cid:15) E F − T , σ L = J F − σ E F − T . (2.21) We recall conservation equations using the Eulerian and Lagrangian for- malisms.
In Eulerian setting, the conservation equations read [36, 31, 37]:
DρDt + ρ div v = 0 , ρ D v Dt − div σ = f E , (cid:15) P : σ + L E = 0 . (2.22 a-c) ρ DU E Dt + div q E = σ : L + w E . (2.23)9quations (2.22 a-c) are balance equations of the mass, linear and angular mo- mentum and (2.23) is the balance equation of the internal energy. The quantity ρ is the mass density (kg m − ), σ is the Cauchy stress (N/m ), f E is the volume force (N/m ), (cid:15) P is the third order permutation tensor also known as the Levi Civita tensor such that ( (cid:15) p σ ) i = ( (cid:15) p ) ijk σ jk , L E is the torque, U E is the den- sity of internal energy, q E is the heat flux density, L := grad v is the gradient deformation and w E is the electromagnetic source term for the heat problem. The electromagnetic force, torque, internal energy and source term are givenby: I f E = ρ E e + j × b + ( grad e ) T p + ( grad b ) T m + ∂ t ( p × b ) + div [ v ⊗ ( p × b )] , (2.24) L E = p × e + ( m − v × p ) × b , L E = p × e + ( m − v × p ) × b . (2.25) ρ DU E Dt = ρ c p ∂ϑ E ∂t , L E = p × e + ( m + v × p ) × b , (2.26) w E = j eff · e eff − m eff · ∂ b ∂t + ρ ∂∂t (cid:18) p ρ (cid:19) · e eff (2.27)where c p and ϑ E are the heat capacity and the temperature, respectively. Equa- tions (2.22 c) and (2.25) suggest that Cauchy stress σ may not be symmetric in presence of electromagnetic fields. However, symmetry may be kept for isotropic materials. Equations (2.22)–(2.23) must be completed by constitutive laws derived fromthe
Clausius-Duhem inequality . In this paper, we assume the following nonlinearconstitutive law for the thermal problem [37]: q E = Q ( ϑ E , grad ϑ E ) = − κ E ( ϑ E ) grad ϑ E (2.28)where κ E is the thermal conductivity tensor (W/mK). The thermo-mechanicalconstitutive law involves the definition of an appropriate objective rate σ ∇ (e.g., Jaumann rate, Truesdell rate or Green-Naghdi rate) which is related tothe material derivative of the Cauchy stress ˙ σ [37] as: σ ∇ = ˙ σ − ψ ( L , F ) = Σ ( σ , L , F , ϑ E , Z E ( τ ≤ t )) , L = ˙ F F − (2.29 a-b)10here the velocity gradient L is related to the rate of the deformation gradient ˙ F through (2.29 b) and Z E ( τ ≤ t ) is the set of internal variables that account for the history of the loadings. In Lagrangian setting, the conservation equations read [36, 31, 37]: ρ − ρJ = 0 , ρ D u Dt − Div
F S = f L , (cid:15) P : ( F SF T ) + L L = 0 . (2.30 a-c) ρ DU L Dt + div q L = − P T : Grad ( F V ) + w L . (2.31)Equations (2.30 a-c) are balance equations of the mass, linear and angular mo-mentum and (2.31) is the balance equation of the internal energy, expressedon the undeformed configuration Ω Mec0 . The quantity ρ is the mass density(kg/m ), S is the second Piola–Kirchhoff stress (N/m ) with S = F − P where P is the first Piola–Kirchhoff stress or nominal stress tensor with P = J σ F − T .The quantity f L is the volume force (N/m ), L L is the torque, U L is the densityof internal energy, q L is the heat flux density and w L is the source term for theheat equation. The electromagnetic force which is a three differential form istransformed as: f L = J f E = ρ L F − T E + J − (( F J ) × ( F B ))+ J (cid:16) F − T Grad ( F − T E ) (cid:17) T ( J − F P )+ (cid:16) F − T Grad ( J − F B ) (cid:17) T ( F − T M )+ J ∂∂t (cid:2) J − ( F P ) × ( F B ) (cid:3) + Div (cid:2) J − V ⊗ [( F P ) × ( F B )] (cid:3) . (2.32)The torque and the internal energy are transformed as: L L = J L E = ( F P ) × (cid:16) F − T E (cid:17) + (cid:16) F − T M eff (cid:17) × ( F B ) , (2.33a) ρ DU L Dt = ρ c p ∂ϑ L ∂t , L E = p × e + ( m + v × p ) × b , (2.33b)11here ϑ L is the temperature expressed on the undeformed configuration andthe source term is obtained using the transformation: w L = J w E = ( F J ) · ( F − T E eff ) (2.34a) − J (cid:16) F − T M eff (cid:17) × (cid:0) ∂ t ( J − F B ) + Grad (cid:0) J − F B (cid:1) V (cid:1) (2.34b) − ρ (cid:18) ∂∂t (cid:18) F P ρ (cid:19) + Grad (cid:18) F P ρ (cid:19) V (cid:19) · (cid:16) F − T E eff (cid:17) . (2.34c) Equations (2.30)–(2.31) must be completed by constitutive laws which relatethe stress tensor to its associated conjugate strain tensor. In this paper, we usethe second Piola–Kirchhoff stress and the Green–Lagrange strain tensors. Weassume the following nonlinear constitutive laws for the thermal and thermo-mechanical problems [37]: q L = J F − Q ( Grad ( ϑ E ◦ ϕ t − ))= J F − κ E F − T (cid:124) (cid:123)(cid:122) (cid:125) κ L Grad ϑ L = − κ L ( ϑ L ) Grad ϑ L , (2.35) S = S VEP ( E , ˙ E , ϑ L , Z L ( τ ≤ t )) , E = 12 ( F T F − ) , (2.36 a-b)where κ L is the thermal conductivity tensor (W/mK) and S VEP is a mapping that represents the visco-elastoplastic constitutive law. Viscosity is reflected through the dependence of the stress on the rate of the Green–Lagrange ten- sor ˙ E , plasticity is accounted for using a set of internal variables Z L ( τ ≤ t ) and the thermo-mechanical aspect is accounted for by the dependency on the temperature ϑ L .
3. Formulations, discretization and linearization
In this section, we derive the simplified multiphysics problem from the fully coupled problem defined in section 2. Using this simplified problem, strong and weak formulations of the multiphysics problems are derived using poten- tials. The weak formulations are then discretized in space using the continuous nally, the resulting nonlinear system of algebraic equations is linearized.
We make the following magnetoquasistatic (MQS) assumptions δ i (cid:39) L sys , i , λ (cid:29) L sys , i . (3.1)In (3.1), δ i := 1 / (cid:112) πf σ E,i µ E,i is the skin depth in the spatial direction i = x, y and z , f is the frequency of the source term, σ E,i and µ E,i are eigenvalues of σ E and µ E , λ is the wavelength corresponding to the frequency f and L sys , i is the characteristic length of the structure in a spatial direction i [38]. The first term of (3.1) explains the presence of eddy currents while the second explains the neglect of electromagnetic waves. Further in this section, we will relax the first condition in (3.1) to δ i = α L sys , i with α which is big thus allowing to neglect the reaction field. Using the MQS assumption, the following MQS problem in the deformedconfiguration can be defined: curl h = j , curl e = − ∂ t b , div b = 0 , (3.2 a-c)together with the constitutive laws: h = ν E ( b ) b = H ( b ) , j = j eff = σ E e eff + j s = σ E ( e + v × b ) + j s . (3.3 a-b) In Lagrangian setting, the MQS problem is governed by Maxwell’s equations[20, 35]:
Curl H eff = J , Curl E eff = − ∂ t B , Div B = 0 (3.4 a-c)completed by the following constitutive laws: H eff = ν L B , J = σ L E eff + J s . (3.5 a-b) .1.2. Conservation equations In addition to the MQS assumption, we assume quasistatic mechanical prob-lem thus neglecting the inertia term in the balance of linear momentum and isotropic magnetic materials therefore restoring the symmetry of the Cauchyand the second Piola–Kirchhoff stress as L E = L L = in (2.22 c) and (2.30 c).Additionally, we neglect mechanical losses in the heat equation. Under theseassumptions, balance equations in the deformed configuration become: DρDt + ρ div v = 0 , div σ + f E = 0 , σ = σ T . (3.6 a-c) ρ c p ∂ϑ E ∂t + div q E = w E (3.7)and the electromagnetic force, torque and electromagnetic losses in (2.24)–(2.27)become: f E = j × b +( grad b ) T m , L E = m × b , w E = j eff · e eff − m eff · ∂ b ∂t . (3.8 a-c) In Lagrangian setting, the conservation equations become: ρ − ρJ = 0 , Div
F S + f L = , S = S T , (3.9 a-c) ρ c p ∂ϑ L ∂t + Div q L = w L . (3.10)where the electromagnetic force and torque are given by: f L = J − (( F J ) × ( F B )) + J (cid:16) F − T Grad ( J − F B ) (cid:17) T ( F − T M ) , (3.11) L L = (cid:16) F − T M eff (cid:17) × ( F B ) , (3.12)and the source term is given by: w L = ( F J ) · ( F − T E eff ) − J ( F − T M eff ) × (cid:2) ∂ t ( J − F B ) + Grad (cid:0) J − F B (cid:1) V (cid:3) . (3.13)14urther in the paper, we consider elasto-plastic materials governed by the fol-lowing constitutive law: S = S EP ( E , ϑ L , Z L ( τ ≤ t )) , E = 12 ( F T F − ) . (3.14 a-b)In this case, the second Piola–Kirchhoff stress does not depend on the rate of Green–Lagrange strain as viscosity is not considered.
The strong formulation of the MQS problem is derived using the so-called magnetic induction conforming formulations [39]. In the deformed configura-tion, the derivation is achieved by verifying equations (3.2 b-c) in the strongsense: b = curl a (cid:39) b s = curl a s , e = − ∂ t a − grad φ (cid:39) − ∂ t a s − grad φ (3.15)where a is the vector potential unknown field used for the eddy currents problem, a s is the source vector potential which can be pre-computed based on the electric current source j s imposed in inductors such as in coils and φ is the unknown scalar potential. The approximation a (cid:39) a s implies the neglect of the reaction field and is valid under the assumption δ i = α L sys , i with α which is big [40]. The electric potential φ is governed by the following problem where (3.16)is obtained by applying the divergence div operator to (3.2 a) and (3.18) isderived from (3.3 b): div j = 0 in Ω Elet , (3.16) j = σ E ( e + v × b ) + j s (3.17) = σ E ( − ∂ t a s − grad φ + v × curl a s ) + j s in Ω Elet , (3.18) φ ( x , t ) = φ D ( x , t ) on Γ Diri , Ele t , (3.19) n E · j = 0 on Γ Neu , Ele t . (3.20)(3.19) and (3.20) are the Dirichlet and the Neumann boundary conditions.
15n the undeformed configuration, the derivation is carried out by verifyingequations (3.4 b-c) in the strong sense: B = Curl A (cid:39) B s = Curl A s = Curl (cid:16) F T a s (cid:17) , E eff = − ∂ t A − Grad Φ (cid:39) − ∂ t A s − Grad
Φ = − F T ∂ t a s − Grad Φ . (3.21)In (3.21), the 1 differential form A s is transformed as A s = F T a s . The electricpotential Φ is therefore governed by the following problem where (3.22) is ob-tained by applying the divergence operator Div to (3.4 a) and (3.24) is derivedfrom (3.5 b):Div J = 0 , in Ω Ele0 , (3.22) J = σ L E eff + J s (3.23) = − σ L ( F T ∂ t a s + Grad
Φ) + J s , in Ω Ele0 (3.24) Φ( x , t ) = Φ D ( x , t ) on Γ Diri , Ele0 , (3.25) n L · J = 0 on Γ Neu , Ele0 . (3.26)where the conductivity tensor is transformed as σ L = J F − σ E ( ϑ E ◦ ϕ − t ) F − T thanks to (2.21). In the deformed configuration, the evolution of the temperature is governedby the following problem derived from (3.7) and (3.8 c) and (2.35): ρc p ∂ϑ E ∂t + div q E = w E in Ω Thet , (3.27) q E = − κ E ( ϑ E ) grad ϑ E in Ω Thet , (3.28) ϑ E ( x ,
0) = ϑ E, ( x ) in Ω The0 , (3.29) ϑ E ( x , t ) = ϑ E,D ( x , t ) on Γ Diri , The t , (3.30) n E · q E = h E ( t )( ϑ E − ϑ E,B ) on Γ conv , The t , (3.31) n E · q E = (cid:15) RE σ RE ( ϑ E − ϑ E,R ) on Γ rad , The t . (3.32)16he source term in (3.27) can be expressed in terms of the potentials as: w E = ( σ E ( ∂ t a s + grad φ − v × curl a s )) (cid:124) (cid:123)(cid:122) (cid:125) w eddy E − j s · ( σ E ( ∂ t a s + grad φ − v × curl a s )) (cid:124) (cid:123)(cid:122) (cid:125) w eddy E − m eff · curl (cid:18) ∂ a s ∂t (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) w hyst E (3.33)where w eddy E and w hyst E represent eddy current and hysteretic losses. Equations (3.29), (3.30), (3.31) and (3.32) represent the initial condition, Dirichlet, con- vective and radiative boundary conditions, respectively. The Lagrangian strong form of the heat problem is given by: ρ c p ∂ϑ L ∂t + Div q L = w L in Ω The0 , (3.34) q L = − κ L ( ϑ L ) Grad ϑ L in Ω The0 , (3.35) ϑ L ( x ,
0) = ϑ L, ( x ) in Ω The0 , (3.36) ϑ L ( x , t ) = ϑ L,D ( x , t ) on Γ Diri , The0 , (3.37) n L · q L = h L ( t )( ϑ L − ϑ L,B ) on Γ conv , The0 , (3.38) n L · q L = (cid:15) RL σ RL ( ϑ L − ϑ L,R ) on Γ rad , The t (3.39)where the thermal conductivity is transformed according to (2.21) as κ L = J F − κ E F − T and the convective heat coefficient h L ( t ) is transformed using Nanson’s formula as h L ( t ) = J | F − T n L | h E ( t ) . The source term in (3.34) can be expressed in terms of the potentials as: w L = ( F − T ( σ L ( ∂ t A s + Grad Φ) − J s (cid:124) (cid:123)(cid:122) (cid:125) − J )) · ( F − T ( ∂ t A s + Grad Φ) (cid:124) (cid:123)(cid:122) (cid:125) − E eff )+ J F − T µ − ( − ν L ) Curl A s (cid:124) (cid:123)(cid:122) (cid:125) M eff (cid:20) ∂ t ( J − F Curl A s ) + Grad (cid:0) J − F Curl A s (cid:1)(cid:21) . (3.40)Equations (3.36), (3.37), (3.38) and (3.39) represent the initial condition, Dirich- let, convective and radiative boundary conditions, respectively. .2.3. The mechanical problem Mechanical fields in the undeformed configuration are governed by the fol-lowing problem:Div ( F S ) + f L = 0 in Ω Mec0 , (3.41) S = S EP ( E , ϑ L , Z L ( τ ≤ t )) in Ω Mec0 , (3.42) E = F T F − in Ω Mec0 , (3.43) u ( x , t ) = u D ( x , t ) on Γ Mec0 ,D . (3.44) n L · ( F S ) = t L on Γ Mec0 ,N . (3.45)In terms of the potential, the force f L in (3.11) is given by: f L = F L ( u , Φ) = J − F ( σ L ( ∂ t A s + Grad
Φ) + J s ) (cid:124) (cid:123)(cid:122) (cid:125) J × F Curl ( A s ) (cid:124) (cid:123)(cid:122) (cid:125) B + J ( F − T Grad ( J − F Curl ( A s ) (cid:124) (cid:123)(cid:122) (cid:125) B )) T ( F − T µ − χ BL J − F T F Curl A s (cid:124) (cid:123)(cid:122) (cid:125) M ) . (3.46)The term t L in (3.45) represents the surface traction applied on part of theboundary Γ Mec0 ,N . The thermo-mechanical constitutive law (3.42)–(3.43) is de-rived from (3.14 a). In this paper, we use the constitutive law described in [30].As a reminder, the total deformation gradient F and the total deformationgradients in the glassy and rubbery states were given by: F = F tg = F tr = + Grad u , (3.47)where the superscript t denotes the total deformation gradient and the super-scripts r and g were used for the rubbery and the glassy states. The totaldeformation gradients of both phases are decomposed as: F tg = F g F f = F eg F pg F f , F tr = F r F p = F er F p , (3.48)where F eg and F pg are the deformation gradients for the elastic and plastic phases in the glassy state, F f is the frozen deformation gradient that represents the temporary deformation which is stored during high temperature shape fixing F r and F p are the elastic and plastic deformation gradients for the rubbery state. The total Cauchy stress was also given by: σ = z g σ g + (1 − z g ) σ r , (3.49)where the temperature-dependent parameter z g is the ratio of the glassy state.Using the second Piola–Kirchhoff stress and Green-Lagrange strain tensors, thefollowing expression of the stress tensor can be derived: S = J F − σF − T = z g J F − σ g F − T + (1 − z g ) J F − σ r F − T = z g J f F f − (cid:16) J g F g − σ g F g − T (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) S g F f − T InnoInnoInno + (1 − z g ) J p F p − (cid:16) J r F r − σ r F r − T (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) S r F p − T . (3.50)In (3.50), S g and S r are the second Piola–Kirchhoff stress tensors defined onthe intermediate configurations [30] by: S g = ( F pg ) − ( λ g tr( E eg ) + 2 µ g E eg ) ( F pg ) − T , (3.51) S r = λ r tr( E er ) + 2 µ r E er . (3.52)The parameters λ g , λ r , µ g and µ r are Lamé parameters for the glassy and the rubbery states. Determinants of deformation gradients of intermediate configu- rations are defined as J i = det F i where the superscript i refer to configurations of the glassy state ( g , eg and pg ) or of the rubbery state ( er and p ). The contribu- tion to the stress due to thermal expansion have been neglected in (3.51)–(3.52). Details on the equations that govern the evolution of internal variables Z L ( τ ≤ t ) := ( z g , F f , F p , F pg ) and the numerical update of the internal vari- ables can be found in [30]. The weak forms of the electromagnetic problem (3.22)–(3.26), the heat prob-lem (3.34)–(3.39) and the mechanical problem (3.41)–(3.45) read [39, 41, 31, 42]:19or each t ∈ I t , find (Φ × ϑ L × u ) ∈ U × V × W such that (cid:90) Ω Ele0 σ L Grad Φ · Grad Φ (cid:48) dΩ + (cid:90) Ω Ele0 σ L ∂ t A s · Grad Φ (cid:48) dΩ = 0 , (3.53) (cid:90) Ω The0 ρ c p ∂ϑ L ∂t · ϑ (cid:48) L dΩ + (cid:90) Ω The0 κ L Grad ϑ L (cid:124) (cid:123)(cid:122) (cid:125) − q L · Grad ϑ (cid:48) L dΩ − InnoInnoInnoInnoInnoInnoInnoInno (cid:90) Ω The0 ( F ( σ L ( ∂ t A s + Grad Φ) − J s )) · ( F − T ( ∂ t A s + Grad
Φ)) (cid:124) (cid:123)(cid:122) (cid:125) w L · ϑ (cid:48) L dΩ + (cid:90) Γ conv , The0 h L ( t )( ϑ L − ϑ L,B ) (cid:124) (cid:123)(cid:122) (cid:125) n L · q L · ϑ (cid:48) L dΓ + (cid:90) Γ rad , The0 (cid:15) RL σ RL ( ϑ L − ϑ L,R ) (cid:124) (cid:123)(cid:122) (cid:125) n L · q L · ϑ (cid:48) L dΓ = 0 (3.54) (cid:90) Ω Mec0 S EP ( u , ϑ L , Z L ) : δ E dΩ − (cid:90) Ω Mec0 F L ( u , Φ) · u (cid:48) dΩ = (cid:90) Γ N0 t L · u (cid:48) dΓ (3.55)holds for all test functions (cid:16) Φ (cid:48) × ϑ (cid:48) L × u (cid:48) (cid:17) ∈ ( U × V × W ) . The force f L = F L ( u , Φ) is given by (3.46) and the dependence on the displacement is achieved through J and F . The function spaces are defined such that U ⊆ H (Ω Ele0 ) , V ⊆ H (Ω The0 ) , W ⊆ H (Ω Mec ) ≡ (cid:0) H (Ω Mec ) (cid:1) and the source term of the electromagnetic problem A s belongs to a subspace of H ( Curl ; Ω
Ele0 ) . The vir- tual Green–Lagrange strain δ E is related to the virtual displacement u (cid:48) through δ E = F T Grad u (cid:48) . The unknown fields Φ , ϑ L and u in (3.53)–(3.55) belong to infinite dimen-sional functional spaces. For numerical simulation, these fields need to be ap-proximated by finite dimensional spaces Φ( x , t ) ≈ ¯Φ( x , t ) , ϑ L ( x , t ) ≈ ¯ ϑ L ( x , t ) , u ( x , t ) ≈ ¯ u ( x , t ) (3.56)20efined by: ¯Φ( x , t ) = N Ele (cid:88) i =1 ¯Φ i ( t ) N Ele i ( x ) , Grad ¯Φ( x ) = N Ele (cid:88) i =1 ¯Φ i ( t ) Grad N i Ele ( x ) , (3.57) ¯ ϑ L ( x , t ) = N The (cid:88) i =1 ¯ ϑ i ( t ) N The i ( x ) , Grad ¯ ϑ L ( x ) = N The (cid:88) i =1 ¯ ϑ i ( t ) Grad N The i ( x ) , (3.58) ¯ u ( x , t ) = N Mec (cid:88) i =1 ¯ u i ( t ) N Mec i ( x ) , Grad ¯ u ( x ) = N Mec (cid:88) i =1 ¯ u i ( t ) Grad N Mec i ( x ) , (3.59)where N Ele , N The and N Mec are the number of nodes of the electromagnetic, thermal and mechanical domains, ¯Φ i , ¯ ϑ i and ¯ u i = (¯ u i,x , ¯ u i,y , ¯ u i,z ) are degrees of freedom and N Ele i , N The i and N Mec i are shape functions for the electromagnetic, thermal and mechanical problems, respectively. Inserting (3.57)–(3.59) into (3.53)–(3.55) leads to the following discrete sys-tem of equations: K Ele (¯ ϑ , ¯ u ) ¯ Φ + F Ele (¯ ϑ , ¯ u ) = , (3.60) M The D ¯ ϑ Dt + K The (¯ ϑ , ¯ u )¯ ϑ + F The ( ¯ Φ , ¯ ϑ , ¯ u ) = , (3.61) K Mec (¯ ϑ , ¯ u , Z L ) + ¯ F Mec ( ¯ Φ , ¯ u ) = . (3.62)where ¯ Φ , ¯ ϑ L and ¯ u are vectors of degrees of freedom and the matrices in (3.60)–213.62) are given by: K Ele = N Ele (cid:88) e =1 L eT (cid:20)(cid:90) Ω e ( B e ) T σ L (¯ u , ¯ ϑ ) B e dΩ e (cid:21) L eT , (3.63) F Ele = N Ele (cid:88) e =1 L eT (cid:20)(cid:90) Ω e ( B e ) T σ L (¯ u , ¯ ϑ ) A es dΩ e (cid:21) , (3.64) M The = N The (cid:88) e =1 L eT (cid:20)(cid:90) Ω e ( N e ) T ρ c p N e dΩ e (cid:21) L eT , (3.65) K The = N The (cid:88) e =1 L eT (cid:20)(cid:90) Ω e ( B e ) T κ L (¯ u , ¯ ϑ ) B e dΩ e (cid:21) L eT , (3.66) F The = N The (cid:88) e =1 L eT (cid:20)(cid:90) Ω e ( B e ) T ¯ w L (¯ u , ¯ ϑ , ¯ Φ ) dΩ e (cid:21) , (3.67) K Mec = N Mec (cid:88) e =1 L eT (cid:20)(cid:90) Ω e ( B e ) T S EP (¯ u , ¯ ϑ , Z L ) dΩ e (cid:21) , (3.68) F Mec = N Mec (cid:88) e =1 L eT (cid:20)(cid:90) Ω e ( N e ) T F L (¯ u , ¯ Φ )dΩ e (cid:21) (3.69)where L e is the gather matrix, N e is the element shape function matrix, B e is composed of the elements of the gradient of N e [43, 31] and where we neglected the boundary terms in (3.53)–(3.55). The recurrent dependence on the displace- ment field in (3.60)–(3.62) and (3.63)–(3.69) results from the transformations between the deformed and undeformed configurations which involve the defor- mation gradient F and its determinant J , and the electric conductivity σ L (¯ u , ¯ ϑ ) is considered to be temperature-dependent. Equations (3.60)–(3.62) can be written as a system of differential algebraicequations: M The
00 0 0 (cid:124) (cid:123)(cid:122) (cid:125) M DDt ¯ Φ ¯ ϑ ¯ u (cid:124) (cid:123)(cid:122) (cid:125) ¯ v + K Ele (¯ ϑ , ¯ u ) ¯ Φ + F Ele (¯ ϑ , ¯ u ) K The (¯ ϑ , ¯ u )¯ ϑ + F The ( ¯ Φ , ¯ ϑ , ¯ u ) K Mec (¯ ϑ , ¯ u , Z L ) + ¯ F Mec ( ¯ Φ , ¯ u ) (cid:124) (cid:123)(cid:122) (cid:125) f (¯ v , Z L ) = (3.70)or M D ¯ v Dt + f (¯ v , Z L ) = , (3.71)22here M is a singular matrix and f = K Ele (¯ ϑ , ¯ u ) ¯ Φ + F Ele (¯ ϑ , ¯ u ) , f = K The (¯ ϑ , ¯ u )¯ ϑ + F The ( ¯ Φ , ¯ ϑ , ¯ u ) , f = K Mec (¯ ϑ , ¯ u , Z L ) + ¯ F Mec ( ¯ Φ , ¯ u ) . (3.72)Equation (3.71) can be discretized in time using the backward Euler integrator: M ¯ v n +1 − ¯ v n ∆ t + f (¯ v n +1 , Z L ) = , (3.73)where ¯ v n +1 = ¯ v ( t n +1 ) with t n +1 = t +( n +1)∆ t and ∆ t which is the time step.After reordering the terms of (3.73), the following nonlinear equations can bederived: M ¯ v n +1 + ∆ t f (¯ v n +1 , Z L ) − M ¯ v n = G (¯ v n +1 , Z L ) − M ¯ v n = H ( ¯Φ n +1 , ¯ ϑ n +1 , ¯ u n +1 , Z L ) − M ¯ v n = (3.74)with the vector function H = ( H , H , H ) defined such that H i ( ¯Φ n +1 , ¯ ϑ n +1 , ¯ u n +1 , Z L ) := G i (¯ v n +1 , Z L ) for i = 1, 2 or 3. Equation (3.74) can be solved using the Newton–Raphson method. To dothis, an iterative schema is used with the following linearization: G (¯ v n +1 , Z L ) (cid:39) G (¯ v n +1 m , Z L ) + (cid:18) ∂ G ∂ ¯ v n +1 (cid:19) ¯ v n +1 m (cid:0) ¯ v n +1 m +1 − ¯ v n +1 m (cid:1) = b RHS (¯ v n +1 m ) − A (¯ v n +1 m ) ∆¯ v n +1 m +1 = , (3.75)where the index m is used to denote the Newton–Raphson iteration. The termsin (3.75) are given by: b RHS = G (¯ v n +1 m , Z L ) = H ( ¯Φ n +1 m , ¯ ϑ n +1 m , ¯ u n +1 m , Z L )= H ( ¯ Φ n +1 m , ¯ ϑ n +1 m , ¯ u n +1 m ) H ( ¯ Φ n +1 m , ¯ ϑ n +1 m , ¯ u n +1 m ) H ( ¯ Φ n +1 m , ¯ ϑ n +1 m , ¯ u n +1 m , Z L ) . (3.76)23he term H depends on the internal variables Z L ( τ ≤ t ) through the second Piola–Kirchhoff stress S . For each quadrature point, the stress is updated using the return mapping described in [30]. The stiffness matrix is given by: (cid:18) ∂ G ∂ ¯ v n +1 (cid:19) ¯ v n +1 m = A = ∂ H ∂ ¯ Φ n +1 ∂ H ∂ ¯ ϑ n +1 ∂ H ∂ ¯ u n +1 ∂ H ∂ ¯ Φ n +1 ∂ H ∂ ¯ ϑ n +1 ∂ H ∂ ¯ u n +1 ∂ H ∂ ¯ Φ n +1 ∂ H ∂ ¯ ϑ n +1 ∂ H ∂ ¯ u n +1 . (3.77)The terms of the tangent stiffness matrix in (3.77) are given by: (cid:18) ∂ H ∂ ¯ Φ n +1 (cid:19) ¯ v n +1 m = (cid:16) K Ele (cid:17) ¯ Φ n +1 m , ¯ ϑ n +1 m , ¯ u n +1 m , (3.78) (cid:18) ∂ H ∂ ¯ ϑ n +1 (cid:19) ¯ v n +1 m = (cid:32) ∂ K Ele ∂ ¯ ϑ n +1 ¯ Φ n +1 m + ∂ F Ele ∂ ¯ ϑ n +1 (cid:33) ¯ Φ n +1 m , ¯ ϑ n +1 m , ¯ u n +1 m , (3.79) (cid:18) ∂ H ∂ ¯ u n +1 (cid:19) ¯ v n +1 m = (cid:32) ∂ K Ele ∂ ¯ u n +1 ¯ Φ n +1 m + ∂ F Ele ∂ ¯ u n +1 (cid:33) ¯ Φ n +1 m , ¯ ϑ n +1 m , ¯ u n +1 m , (3.80) (cid:18) ∂ H ∂ ¯ Φ n +1 (cid:19) ¯ v n +1 m = ∆ t (cid:32) ∂ F The ∂ ¯ Φ n +1 (cid:33) ¯ Φ n +1 m , ¯ ϑ n +1 m , ¯ u n +1 m , (3.81) (cid:18) ∂ H ∂ ¯ ϑ n +1 (cid:19) ¯ v n +1 m = (cid:32) M The + ∆ t (cid:32) K The + ∂ K The ∂ ¯ ϑ n +1 ¯ ϑ n +1 m + ∂ F The ∂ ¯ ϑ n +1 (cid:33)(cid:33) , (3.82) (cid:18) ∂ H ∂ ¯ u n +1 (cid:19) ¯ v n +1 m = ∆ t (cid:32) ∂ K The ∂ ¯ u n +1 ¯ ϑ n +1 m + ∂ F The ∂ ¯ u n +1 (cid:33) ¯ Φ n +1 m , ¯ ϑ n +1 m , ¯ u n +1 m , (3.83) (cid:18) ∂ H ∂ ¯ Φ n +1 (cid:19) ¯ v n +1 m = (cid:32) ∂ F Mec ∂ ¯ Φ n +1 (cid:33) ¯ Φ n +1 m , ¯ u n +1 m , (3.84) (cid:18) ∂ H ∂ ¯ ϑ n +1 (cid:19) ¯ v n +1 m = (cid:32) ∂ K Mec ∂ ¯ ϑ n +1 (cid:33) ¯ ϑ n +1 m , ¯ u n +1 m , (3.85) (cid:18) ∂ H ∂ ¯ u n +1 (cid:19) ¯ v n +1 m = (cid:32) ∂ K Mec ∂ ¯ u n +1 + ∂ F Mec ∂ ¯ u n +1 (cid:33) ¯ Φ n +1 m , ¯ ϑ n +1 m , ¯ u n +1 m . (3.86)In (3.86), the first term is given by (cid:32) ∂ K Mec ∂ ¯ u n +1 (cid:33) ¯ v n +1 m = N Mec (cid:88) e =1 L eT (cid:20)(cid:90) Ω e ( B e ) T (cid:18) ∂ S EP ∂ E : ∂ E ∂ ¯ u n +1 (cid:19) dΩ e (cid:21) . (3.87)24he Jacobian of the mechanical problem ∂ S EP /∂ E in (3.86) is also updated using the return mapping algorithm described in [30]. The pseudocode in Algo- rithm 1 illustrates the flow of the numerical code used to solve the multiphysics problem. Algorithm 1
Pseudocode of the multiphysics SMP problem
INPUT:
Mesh, current source I s ( t ) and constitutive laws. OUTPUT:
Mechanical, thermal and electromagnetic fields, and Joule losses procedure
Multiphysics problem t ← t , initialization of the physica fields T | t = T and u | t = u , for ( n ← To N TS ) do (cid:46) the time loop ( index n ) Innocent Ii for ( m ← To N MNR ) do (cid:46) the NR loop of the overall problem for ( l ← To N RM ) do (cid:46) the NR loop of the return mapping Pass internal variables as input,Update the mechanical law S and ∂ S EP /∂ E using the RM end for Update and compute the Jacobian,Assemble (cid:0) ∂ G /∂ ¯ v n +1 (cid:1) ¯ v n +1 m from (3.77) and G (¯ v n +1 m ) from (3.76),Solve the system (cid:0) ∂ G /∂ ¯ v n +1 (cid:1) ¯ v n +1 m ∆ v n +1 m +1 = − G (¯ v n +1 m , Z L ) ,compute the residual r ( v n +1 m ) = G (¯ v n +1 m ) , if ( || r ( v n +1 m ) || ≤ ε tol ) then Exit the nonlinear loop of the overall problem else
Do another NR iteration for the overall problem end ifend for
Go to the next time step t ← t + ∆ t end forend procedure . Numerical Tests This section is devoted to the numerical testing of the electro-thermo-mechanical problem. A set of numerical tests similar to the ones developed in [30] are herein proposed. Whereas the authors in [30] considered ideal and non-ideal shape memory polymer materials, the main focus of this paper is on the fully coupled problem. Therefore, in Section 4.1 we consider an ideal and a non-ideal shape memory polymer single element similar to the one in [30] for the validation of the thermomechanical problem. In Section 4.2 we only consider an ideal shape memory polymer stent for the fully coupled electro-thermo-mechanical problem.
The two considered tests are: • the uniaxial tests on a × × mm single-element cube (SEC), • the simulation of a cylindrical vascular stent (CVS) similar to the one de- scribed in [30]. The stent has the same dimensions and material properties but without the small holes. Symbol Value Unit E r E g
771 MPa ν r ν g R pg
10 MPa h ∆ θ
30 (SEC) – 5 (CVS) K θ t
350 (SEC) – 344 (CVS) K w c (ideal) (ideal) 1 (ideal) (ideal) – (ideal) c p Table 1: Model parameters of the mechanical problem Material properties listed in Table 1 are used for both cases. The software
GetDP [44] was used to solve the fully coupled problem based on a total La- grangian formulation.
Results of the thermo-mechanical model developed in [30] are reproduced.
This model consisted of a temperature-dependent elasto-plastic model with a constant temperature field imposed for all Gauss points at any given time instant t . Results of the single element are reported in Figure 2 for a high-temperature fixing similar to the one used in [30]. The material is progressively deformed at material is then unloaded at 200K before re-heating it up to 400K to trigger the shape-recovery (see Test 1 in Fig. 6 of [30]). The results reported in Figure 2 conform to those obtained in [30].
The results of the coupled problem are presented below. A description of the mechanical, thermal and electromagnetic problems is followed by the pre- sentation of numerical results of the shape memory polymer. A best design can be obtained by choosing material properties for the thermal problem (thermal conductivity, mass density and heat capacity) that maintain a homogeneous temperature field in the stent, in order to avoid the appearance of regions with different phases during the recovery step. A non dimensionalization analysis of the thermal problem carried out in Section 4.2.2 facilitates this design. How- ever, the control of the temperature is complicated by the dependence of the mechanical stress on the temperature-dependent ratio of the glassy state z g ( θ L ) . The geometry of Figure 3 is used for the mechanical problem of the cylin- drical vascular stent. . − . ε yy (%) T (K) σ yy ( M P a ) . . Time (s) ε yy ( % ) Time (s) σ yy ( M P a ) . . ε yy (%) σ yy ( M P a )
200 300 40005
Temperature (s) σ yy ( M P a )
200 300 40000 . . T (K) ε yy ( % ) Figure 2: Results of the single element. Top-left: temperature–strain–stress curve. Top-right:strain versus time curve. Middle-left: stress versus time curve. Middle-right: strain versusstress curve. Bottom-left: temperature versus stress curve. Bottom-right: temperature versusstrain curve. Blue curves correspond to the ideal case while red curves correspond to thenon-ideal case. Top D u Bottom D Figure 3: Dirichlet boundary conditions for the mechanical problem. Zero displacement u Bottom D ( t ) = 0 is imposed on the bottom section and a time-dependent displacement u Top D ( t ) similar to the one used for Test 1 of Fig. 6 in [30] is imposed on the top line. We simulate the insertion of a vascular shape memory polymer stent in a vein of the arm. The stent contains electric particles that can react to electromag- netic source fields produced by a coil wrapped around the arm by producing heat by the Joule effect. For the sake of simplicity, we consider the resulting shape memory polymer composite to be homogeneous with homogenized macroscopic material properties, thus ignoring the multiscale nature of the composite.
The mechanical problem is similar to the one in [30], with the temperature field obtained by solving the thermal problem with the source generated by the eddy current losses. In the following, we define the electromagnetic and the thermal problems.
The temperature can be controlled by an electromagnetic field generated by29 coil crossed by a current denoted I s ( t ) . For all problems studied herein, weconsider a single frequency source I s ( t ) = I ( t )( a + b sin( ω t )) = I ( t )( a + b sin(2 π f t )) , (4.1)where I ( t ) (A) is piecewise, linear, time-dependent amplitude of the electriccurrent, ω is the angular velocity and f the frequency of the signal. The designparameters for the electromagnetic and thermal problems are the amplitude ofthe current, the frequency and the material properties: the electric conductivity σ (S/m), the magnetic permeability µ (H/m), the mass density ρ , the heatcapacity c p and the thermal conductivity κ . In all our applications, we considerfrequencies small than 1000Hz, σ = 10 S/m and µ = µ µ rel with µ rel = 20 ,which corresponds to the wavelength λ and skin depth δ : λ = cf = 1 √ µ(cid:15)f ≈ km , δ = (cid:114) µσω ≈ mm . (4.2)The wavelength is very large compared to the dimensions of the structure (typically 20mm for the length and 1mm for the thickness) that the quasistatic assumption can be made [45]. Likewise, the skin depth is large compared to the dimensions of the stent that the eddy currents resulting from the reaction field can be neglected. Figure 4 illustrates the geometry used for the coupled problem. To determine the magnetic induction source b s ( t ) for the electromagneticproblem, we consider a coil with a very large number of turns. The valueof the magnetic field h s ( t ) and the magnetic induction b s ( t ) in the coil arehomogeneous and given by [46] h s ( t ) = h s ( t ) e z = NL I s ( t ) e z , b s ( t ) = b s ( t ) e z = µ h s ( t ) = µ NL I s ( t ) e z (4.3)where N is the number of turns, L the length of the coil, µ the magnetic per-meability of the material and e z the direction oriented along the axis of thecoil. From the Gauss magnetic law div b s = 0 , a source vector potential a s ( t ) can be derived from the magnetic induction b s ( t ) as b s ( t ) = curl a s ( t ) . In thecomputational domain of the stent, a possible vector potential that satisfies this30 igure 4: Geometry and mesh used for the coupled problem. Top: The cylindrical cardiovas-cular stent surrounded by an exciting coil. Middle: Mesh of the stent and the coil. Bottom:Mesh of the stent, the coil and the surrounding air. The enclosing box is used to bound thecomputational domain for the electromagnetic problem assumed to be unbounded. a s ( x, y, t ) = 0 . b s ( t )( − y, x,
0) = 0 . µ NL I s ( t )( − y, x, (4.4)with x = X + u x and y = Y + u y where X and Y are the coordinates expressed in the undeformed configuration and u x and u y are components of the displacement u = ( u x , u y , u z ) . From (4.3) and (4.4), it can be noted that the magnetic field h s ( t ) and the magnetic induction b s ( t ) in the coil do not depend on spatial coordinates whereas the vector potential a s ( t ) depends on spatial coordinates x and y . This vector potential is a one differential form that can be transformed as A s = F T a s thus leading to the source in (2.15). Table 2 contains model parameters of the electromagnetic and thermal prob- lems.
Symbol Value Unit I ( t ) electric current waveform A f σ S/m µ r
20 – ρ
270 kg m − c p
10 kg m K − s − (ideal) k (ideal) (ideal) 237 (ideal) (ideal) W m − K − (ideal) h
500 W m − K − N L Table 2: Model parameters of the electromagnetic and thermal problems
Defining the thermal problem resulting from the deployment of the actual stent is challenging. Though it is easy to control the temperature of the device during the first three stages ( loading at high temperature, cooling followed by nloading/insertion of the stent) most of which are done outside the human body, the last step, the recovery , necessitates controlling the temperature us- ing electromagnetic fields. In this paper, we simulate the control of the entire deployment process using the electromagnetic fields. During the last step, different modes of heat exchange can be considered: (1) heat conduction in the stent, at the interface of the stent and the surrounding tissue and in the tissue itself and (2) forced convection at part of the boundary of surface of the stent in contact with the blood flowing in the vein. The surface of the stent in contact with the tissue/blood varies during the process of recovery and its detection would necessitate consideration of contact mechanics. For the sake of simplicity, we only consider forced convection.
Finally, thanks to the assumptions made of the electromagnetic and thermal problems, all three problems can only be solved on the computational domain of the stent thus neglecting the surrounding environment.
Results of the coupled problem are herein reported. As mentioned earlier, the main difference between this section and section 4.1 lies in the use of a temperature field obtained by solving the heat equation on a moving domain with the source obtained by solving the electromagnetic problem instead of a priori imposing a temperature field at each time instant t . Figures 5 and 6 show the displacement u , the current density j , Joule losses and the temperature T at the instances t = 4 . × − s and t = . × − s. This can play an important role in the design of the stent, especially for the computation of the temperature. Indeed, the high depen- dency of the ratio of the glassy states on the temperature z g ( θ L ) necessitates selecting electromagnetic and thermal loadings as well as thermal material prop- erties that allow for a quick diffusion of the heat sources throughout the stent, to avoid inhomogeneities of temperature that would cause different regions of the stent to be in different phases (rubbery/glassy) during the recovery step. Another issue concerns the use of the Newton–Raphson method to solve the .
000 0 .
007 0 . Displacement (m) .
000 0 .
003 0 . Displacement (m) . e +06 1 . e +06 1 . e +06 Current density (A/m ) . e +06 2 . e +06 2 . e +06 Current density (A/m ) . e +08 1 . e +08 2 . e +08 Joule Losses (W) . e +08 4 . e +08 4 . e +08 Joule Losses (W)
Figure 5: Physical fields during the recovery. The left images correspond to t = 4 . × − sand the right images correspond to t = 4 . × − s. Top : displacement u , middle :current density J , bottom : Joule losses w L . . e +02 3 . e +02 3 . e +02 T (K) . e +02 3 . e +02 3 . e +02 T (K)
Figure 6: Temperature at t = 4 . × − s (left) and t = 4 . × − s (right). nonlinear coupled problem. Considerably large and inhomogeneous increments of temperature computed especially during the first nonlinear iterations of the Newton–Raphson scheme may lead to inhomogeneities of temperature and cause slow convergence in the recovery process.
To avoid inhomogeneities of temperature in the recovery step, we developed the following normalization process, which makes the problem well conditioned.
The process starts with the linearized version of the heat equation (3.34): ρ L c p ∂θ L ∂t + Div X [ κ L Grad X θ L ] = − w L ( φ, θ L , u ) . (4.5)A new coordinate system ( τ , η ) is introduced as: t = T c τ , dt = T c dτ , ∂ ( · ) ∂t = 1 T c ∂ ( · ) ∂τ ,X i = L c η i , dX i = L c dη i , ∂ ( · ) ∂X i = 1 L c ∂ ( · ) ∂η i , θ L = θ c ¯ θ L (4.6)where ¯ θ , T c , L c , τ , and η i are the characteristic temperature, the characteristictime, the characteristic length, the dimensionless temporal and spatial coordi-nates, respectively. The derivatives in (4.5) are transformed as: ∂θ L ∂t = 1 T c ∂ (cid:0) θ c ¯ θ L (cid:1) ∂τ = θ c T c ∂ ¯ θ L ∂τ , ∂θ L ∂X i = 1 L c ∂ (cid:0) θ c ¯ θ L (cid:1) ∂η i = θ c L c ∂ ¯ θ L ∂η i ,∂ θ L ∂X i = 1 L c ∂ (cid:18) ∂θ L ∂x i (cid:19) ∂η i = θ c L c ∂ ¯ θ L ∂η i (4.7)35 · − Displacement u y (mm) F y ( N ) · − Displacement u y (mm) F y ( N )
320 340 360051015 T (K) D i s p l a ce m e n t u y ( mm )
320 340 360051015 T (K) D i s p l a ce m e n t u y ( mm ) Figure 7: Results of the coupled problem. Top : Force versus displacement curves. Com-puted temperature (left) and imposed temperature (right). Bottom: Displacement versustemperature curves. Computed temperature (left) and imposed temperature (right). which leads to the dimensionless heat equation : ∂ ¯ θ L ∂τ + T c ¯ κ L ρ L c p L c Div η (cid:2) Grad η ¯ θ L (cid:3) = − T c ρ L c p L c ¯ w L ( ¯ φ, ¯ θ L , ¯ u ) (4.8)where the thermal conductivity was assumed constant and barred quantities aredefined in the new coordinate system. For the first two terms to be of the sameorder of magnitude, i.e., for the temperature to have enough time to diffuse inthe stent, the material properties must be chosen such that T c = ρ L c p L c ¯ κ L . (4.9)Results of the coupled problem are reported in Figures 7-8 for a stent with the high-temperature fixing and slightly different material properties as those
20 340 36002468 · − T (K) F y ( N )
320 340 36002468 · − T (K) F y ( N ) · − t (s) T ( K ) t (s) T ( K ) Figure 8: Results of the coupled problem. Top : Force versus temperature curves. Computedtemperature (left) and imposed temperature (right). Bottom : Temperature versus timecurves. Computed temperature (left) and imposed temperature (right). reported in [30]. In the case of the imposed temperature, the material is pro- gressively deformed at 350K, then cooled down to 320K while the deformation is maintained constant, then unloaded at 320K, and finally re-heated up to 350K to trigger shape-recovery. We mimick the same trajectory of the temperature by changing material properties, the frequency and the amplitude of the excitation source. Results of the coupled problem are different from the ones obtained with the imposed temperature. An optimal control of the temperature using the source current I s ( t ) and the geometry of the coil as control parameters can allow to prescribe a temperature trajectory convenient for surgical purposes. . Conclusions In this paper, the deployment of a vascular shape memory polymer stent in a vein of an arm is simulated. The temperature field used in the thermo- mechanical model of the stent is controlled by solving for electromagnetic fields generated by a coil wrapped around the arm. The controllability of the tem- perature depends on the choice of the electromagnetic source field determined by the amplitude and the excitation frequency of the current flowing through the coil, and on the material properties used for the thermal problem. An ini- tial design of the stent which allows for the diffusion of heat and leads to a homogeneous distribution of temperature during the recovery step is proposed.
The optimal control of the temperature of the devices can further be carried out, thus allowing the device to follow a prescribed temperature trajectory that might be convenient for surgical purposes.
Acknowledgments
The authors would like to thank Dr. Elisa Boatti at Georgia Institute of
Technology, Prof. Ludovic Noels and Miguel Pareja Mu ˜ noz at the University of Liège for fruitful discussions on mechanical models of shape memory polymers.
They would also like to thank Prof. Christophe Geuzaine at the University of
Liège for discussions about the implementation in GetDP. The first author is particularly indebted to Dr. François Henrotte at the University of Liège for the discussions on electromagnetic formulations under large deformations. During the time the research was carried out, Innocent Niyonzima was a postdoctoral
Fellow with the Belgian American Educational Foundation (BAEF). He is also partially supported by an excellence grant from Wallonie-Bruxelles International (WBI).
References [1] An introduction to biomaterials, edu/research/tutorials/introbiomat.html , accessed: 2017-12-22.
Business Media, 2010. [3] H. R. Rezaie, L. Bakhtiari, A. Öchsner, Biomaterials and their applications,
Springer, 2015. [4] W. M. Huang, B. Yang, Y. Q. Fu, Polyurethane shape memory polymers,
CRC Press, 2011. [5] W. Small, P. R. Buckley, T. S. Wilson, W. J. Benett, J. Hartman, D. Sa- loner, D. J. Maitland, Shape memory polymer stent with expandable foam: a new concept for endovascular embolization of fusiform aneurysms, IEEE
Transactions on Biomedical Engineering 54 (6) (2007) 1157–1160. [6] D. J. Maitland, W. Small, J. M. Ortega, P. R. Buckley, J. Rodriguez,
J. Hartman, T. S. Wilson, Prototype laser-activated shape memory polymer foam device for embolic treatment of aneurysms, Journal of biomedical optics 12 (3) (2007) 030504–030504. [7] C. M. Yakacki, R. Shandas, C. Lanning, B. Rech, A. Eckstein, K. Gall, Un- constrained recovery characterization of shape-memory polymer networks for cardiovascular applications, Biomaterials 28 (14) (2007) 2255–2263. [8] G. M. Baer, T. S. Wilson, W. Small, J. Hartman, W. J. Benett, D. L.
Matthews, D. J. Maitland, Thermomechanical properties, collapse pressure, and expansion of shape memory polymer neurovascular stent prototypes,
Journal of Biomedical Materials Research Part B: Applied Biomaterials
90 (1) (2009) 421–429. [9] S. H. Ajili, N. G. Ebrahimi, M. Soleimani, Polyurethane/polycaprolactane blend with shape memory effect as a proposed material for cardiovascular implants, Acta biomaterialia 5 (5) (2009) 1519–1530. [10] L. Yahia, Shape memory polymers for biomedical applications, Elsevier, drug loading of shape-memory polymer networks–effect on their functional- ities, European Journal of Pharmaceutical Sciences 41 (1) (2010) 136–147. [12] K. Nagahama, Y. Ueda, T. Ouchi, Y. Ohya, Biodegradable shape-memory polymers exhibiting sharp thermal transitions and controlled drug release,
Biomacromolecules 10 (7) (2009) 1789–1794. [13] J. Fish, W. Chen, Modeling and simulation of piezocomposites, Computer methods in applied mechanics and engineering 192 (28-30) (2003) 3211– [14] M. Elhadrouz, T. B. Zineb, E. Patoor, Finite element analysis of a multi- layer piezoelectric actuator taking into account the ferroelectric and ferroe- lastic behaviors, International journal of engineering science 44 (15) (2006) [15] P. I. Anderson, A. J. Moses, H. J. Stanbury, Assessment of the stress sensi- tivity of magnetostriction in grain-oriented silicon steel, IEEE transactions on magnetics 43 (8) (2007) 3467–3476. [16] M. Khalaquzzaman, B. Xu, S. Ricker, R. Müller, Computational homoge- nization of piezoelectric materials using FE to determine configurational forces, Technische Mechanik 32 (1) (2012) 21–37. [17] S. Kuznetsov, J. Fish, Mathematical homogenization theory for electroac- tive continuum, International Journal for Numerical Methods in Engineer- ing 91 (11) (2012) 1199–1226. [18] O. Perevertov, J. Thielsch, R. Schäfer, Effect of applied tensile stress on the hysteresis curve and magnetic domain structure of grain-oriented transverse Fe-3% Si steel, Journal of Magnetism and Magnetic Materials 385 (2015) [19] P. L. Bishay, S. N. Atluri, Computational Piezo-Grains (CPGs) for a highly-efficient micromechanical modeling of heterogeneous piezoelectric– (2015) 311–328. [20] A. C. Eringen, G. A. Maugin, Electrodynamics of continua I: foundations and solid media, Springer Science & Business Media, 2012. [21] Y.-H. Pao, Electromagnetic forces in deformable continua, in: Mechanics today, Vol. 4, 1978, pp. 209–305. [22] R. Ogden, Incremental elastic motions superimposed on a finite deforma- tion in the presence of an electromagnetic field, International Journal of
Non-Linear Mechanics 44 (5) (2009) 570–580. [23] P. Saxena, On the general governing equations of electromagnetic acoustic transducers, Archive of Mechanical Engineering 60 (2) (2013) 231–246. [24] M. Stiemer, J. Unger, B. Svendsen, H. Blum, An arbitrary Lagrangian Eule- rian approach to the three-dimensional simulation of electromagnetic form- ing, Computer Methods in Applied Mechanics and Engineering 198 (17) (2009) 1535–1547. [25] B. E. Abali, A. F. Queiruga, Theory and computation of electromagnetic fields and thermomechanical structure interaction for systems undergoing large deformations, https://arxiv.org/abs/1803.10551 , accessed: 2018- [26] G. Ethiraj, C. Miehe, Multiplicative magneto-elasticity of magnetosensitive polymers incorporating micromechanically-based network kernels, Interna- tional Journal of Engineering Science 102 (2016) 93–119. [27] C. Miehe, D. Vallicotti, S. Teichtmeister, Homogenization and multiscale stability analysis in finite magneto-electro-elasticity. Application to soft matter EE, ME and MEE composites, Computer Methods in Applied Me- chanics and Engineering 300 (2016) 294–346. particulate magnetoactive composites, Journal of Engineering Materials and Technology 140 (1) (2018) 011003. [29] L. Homsi, L. Noels, A discontinuous Galerkin method for non-linear electro- thermo-mechanical problems: application to shape memory composite ma- terials, Meccanica (2017) 1–45. [30] E. Boatti, G. Scalet, F. Auricchio, A three-dimensional finite-strain phe- nomenological model for shape-memory polymers: Formulation, numerical simulations, and comparison with experimental data, International Journal of Plasticity 83 (2016) 153–177. [31] P. Wriggers, Nonlinear finite element methods, Springer Science & Business
Media, 2008. [32] P. L. Penfield Jr, L. Chu, H. A. Haus, Electrodynamics of moving media,
Tech. rep., Research Laboratory of Electronics (RLE) at the Massachusetts
Institute of Technology (MIT) (1963). [33] R. M. Fano, L. J. Chu, R. B. Adler, Electromagnetic fields, energy, and forces, MIT Press, 1968. [34] J. D. Jackson, Classical Electrodynamics, 3rd Edition, Wiley, 1998. [35] A. B. de Castro, D. Gómez, P. Salgado, Mathematical models and numer- ical simulation in electromagnetism, Vol. 74, Springer, 2014. [36] J. C. Simo, T. J. Hughes, Computational inelasticity, Vol. 7, Springer Sci- ence & Business Media, 2006. [37] T. Belytschko, W. K. Liu, B. Moran, K. Elkhodary, Nonlinear finite ele- ments for continua and structures, John Wiley & Sons, 2013. [38] R. Hiptmair, O. Sterz, Current and voltage excitations for the eddy current model, International Journal of Numerical Modelling: Electronic Networks,
Devices and Fields 18 (1) (2005) 1–21.
Complementarity, Edge Elements, Academic Press, 1998. [40] R. Scorretti, R. V. Sabariego, L. Morel, C. Geuzaine, N. Burais, L. Nicolas,
Computation of induced fields into the human body by dual finite element formulations, IEEE Transactions on Magnetics 48 (2) (2012) 783–786. [41] F. Bachinger, U. Langer, J. Schöberl, Numerical analysis of nonlinear mul- tiharmonic eddy current problems, Numerische Mathematik 100 (4) (2005) [42] T. J. Hughes, The finite element method: linear static and dynamic finite element analysis, Courier Corporation, 2012. [43] J. Fish, T. Belytschko, A first course in finite elements, John Wiley, 2007. [44] P. Dular, C. Geuzaine, F. Henrotte, W. Legros, A general environment for the treatment of discrete problems and its application to the finite element method, IEEE Transactions on Magnetics 34 (5) (1998) 3395–3398. [45] A. A. Rodríguez, A. Valli, Eddy Current Approximation of Maxwell Equa- tions: Theory, Algorithms and Applications, Vol. 4, Springer Science &
Business Media, 2010. [46] The electromagnetic fields generated by a long solenoid, http:// hyperphysics.phy-astr.gsu.edu/hbase/magnetic/solenoid.html , ac- cessed: 2018-06-04.517