High-Performance Numerical Simulation of Biodegradation Process with - PowerPoint PPT Presentation

High-Performance Numerical Simulation of Biodegradation Process with Moving Boundaries FreeFEM Days, 11th Edition Mojtaba Barzegari, Liesbet Geris Biomechanics Section, Department of Mechanical Engineering, KU Leuven December 2019

0 Our Research Group ◮ Supervisor: Prof. Ir. Liesbet Geris ◮ Research profile: Computational Tissue Engineering, Computational Biomechanics, Computational Biology, Computational Genomics 1

1 Outline 1 Introduction 2 Mathematical Model 3 Computational Model and Parallelization 4 Simulation Results 5 Performance Analysis 2

1 Reaction-Diffusion Systems with Moving Boundaries ◮ Stefan problems ◮ Diffusion-controlled interface ◮ Diffusion and reaction lead to the change of domain geometry ◮ Degradation is an example of such a system 3

1 Biodegradation Process ◮ Dissolution of the bulk material ◮ Formation of a protective film ◮ Effect of ions in the medium 4

1 A Sample Application ◮ Hip joint replacement implants ◮ Tuning the degradation parameters to the rate of bone growth 5

2 Outline 1 Introduction 2 Mathematical Model 3 Computational Model and Parallelization 4 Simulation Results 5 Performance Analysis 6

2 Chemistry of Biodegradation Some of the chemical reactions: Mg → Mg 2+ + 2e − 2H 2 O + 2e − → H 2 + 2OH − Mg 2+ + 2OH − k 1 − → Mg(OH) 2 → Mg 2+ + 2Cl − + 2OH − Mg(OH) 2 + 2Cl − k 2 − 7

2 Reaction-Diffusion Equations x ∈ Ω ⊂ R 3 C Mg = C Mg ( x, t ) , C Film = C Film ( x, t ) � � ∂C Mg C Film � � D e + k 2 C Film [Cl] 2 = ∇ • Mg • ∇ C Mg − k 1 C Mg 1 − ∂t [Film] max � � ∂C Film C Film − k 2 C Film [Cl] 2 = k 1 C Mg 1 − ∂t [Film] max �� � � C Film C Film ǫ D e Mg = D Mg 1 − + [Film] max [Film] max τ 8

2 Level Set Method x ∈ Ω ⊂ R 3 Implicit signed distance function φ = φ ( x, t ) − → ∂φ ∂t + V • ∇ φ + v |∇ φ | = bκ |∇ φ | � �� � � �� � � �� � External velocity field Normal direction motion Curvature - dependent term 9

2 Coupling Mass Transfer and Level Set ∂φ ∂t + v |∇ φ | = 0 Rankine-Hugoniot: { J ( x, t ) − ( c sol − c sat ) v( x, t ) } · n = 0 D e Mg ∇ n C Mg − ([Mg] sol − [Mg] sat ) v = 0 D e Mg ∇ n C Mg ∂φ ∂t − |∇ φ | = 0 [Mg] sol − [Mg] sat 10

3 Outline 1 Introduction 2 Mathematical Model 3 Computational Model and Parallelization 4 Simulation Results 5 Performance Analysis 11

3 Weak Formulation Rewriting the diffusion-reaction PDE: ∂u ∂t = ∇ • ( D • ∇ u ) − k 1 bu + k 2 pq 2 Defining trial and test function space: � � u ( x , t ) | x ∈ Ω , t > 0 , u ( x , t ) ∈ H 1 (Ω) , and ∂u S t = ∂n = 0 on Γ � � v ( x ) | x ∈ Ω , v ( x ) ∈ H 1 (Ω) , and v ( x ) = 0 on Γ V = 12

3 Weak Formulation cont. ∂u ∂t v = ∇ • ( D • ∇ u ) v − k 1 buv + k 2 pq 2 v ∀ v ∈ V Integrate over the whole domain: � ∂u � � � Ω k 2 pq 2 vdω ∂t vdω = Ω ∇ • ( D • ∇ u ) vdω − Ω k 1 buvdω + Ω Integration by part, Green’s divergence theory, Backward Euler scheme: u − u n � � � � � Γ Dv • ∂u Ω k 2 pq 2 vdω vdω = ∂ndγ − Ω D • ∇ u • ∇ vdω − Ω k 1 buvdω + ∆ t Ω 13

3 Weak Formulation cont. � � � � � Ω ∆ tk 2 pq 2 vdω Ω u n vdω + Ω uvdω + Ω ∆ tD • ∇ • u ∇ vdω + Ω ∆ tk 1 buvdω = � By defining a linear functional ( f, v ) = Ω fvdω ( u, v )[1 + ∆ tk 1 b ] + ∆ t ( D ∇ u, ∇ v ) = ( u n , v ) + ∆ t ( f n , v ) 1 multiplying to a new coefficient α = 1+∆ tk 1 b ( u, v ) + α ∆ t ( D ∇ u, ∇ v ) = α ( u n , v ) + α ∆ t ( f n , v ) 14

3 Discretization Scheme � � V h = span { ψ i } i ∈I s I s = { 0 , . . . , N } Using 1st order Lagrange polynomials as basis functions N N u n = � � c n u = c j ψ j ( x ) , j ψ j ( x ) j =0 j =0 N N N � � � ( ψ i , ψ j ) c n j + ∆ t ( f n , ψ i ) ( ψ i , ψ j ) c j + α ∆ t ( ∇ ψ i , D ∇ ψ j ) c j = j =0 j =0 j =0 15

3 Discretization Scheme cont. A linear system of equations � A i,j c j = b i j A i,j = ( ψ i , ψ j ) + α ∆ t ( ∇ ψ i , D ∇ ψ j ) N � α ( ψ i , ψ j ) c n j + α ∆ t ( f n , ψ i ) b i = j =0 16

3 Discretization Scheme cont. Final form as implemented in FreeFEM ( M + α ∆ tK ) c = αMc 1 + α ∆ tf M = { M i,j } , M i,j = ( ψ i , ψ j ) , i, j ∈ I s K = { K i,j } , K i,j = ( ∇ ψ i , D ∇ ψ j ) , i, j ∈ I s f = { f i } , f i = ( f ( x , t n ) , ψ i ) , i ∈ I s c = { c i } , i ∈ I s c 1 = { c n i } , i ∈ I s 17

3 Level Set Implementation ◮ Penalization for interface BCs ◮ Computing ∇ n C Mg correctly ◮ Problem of oscillation ◮ Too flat or too steep gradients ◮ Nightmare of re-distancing (P. Bajger et al. 2017) 18

3 Computational Mesh ◮ Eulerian mesh ◮ Generated using Netgen in SALOME platform ◮ Adaptively refined on the moving interface 19

3 Parallelization ◮ Message Passing Interface ◮ Distributed numerical integration (assigning a number in the range of [0, MPI Size-1] to each element) ◮ MUMPS multifrontal direct solver

4 Outline 1 Introduction 2 Mathematical Model 3 Computational Model and Parallelization 4 Simulation Results 5 Performance Analysis 21

4 Release of Ions and Degradation - Simple Screw 22

4 Film Formation - Simple Screw Formation of the protective film on the interface of material-medium 23

4 Release of Ions and Degradation - Porous Structure Trimmed view of the computational mesh Formation of the protective film 24

4 Quantitative Results Measuring mass loss: ◮ Direct weight reduction ◮ Side products evolution Using level set output for calculating mass loss � � Mg lost = Ω + ( t ) Mg solid d V − Ω + (0) Mg solid d V 0 Ω + ( t ) = { x : φ ( x , t ) ≥ 0 } Simulation and experimental setup 25

4 Mass Loss and Evolving Side Products Film formation and the comparison of predicted and experimental mass loss, measured by the evolved hydrogen gas 26

5 Outline 1 Introduction 2 Mathematical Model 3 Computational Model and Parallelization 4 Simulation Results 5 Performance Analysis 27

5 Problem Size ◮ Same setup as the model for calibration and validation ◮ DOF: 144k ◮ Elements: 831k (P1) 28

5 Domain Decomposition Two different approaches for domain decomposition. Execution time per time step Colors show different mesh regions assigned to different MPI cores. 29

5 Weak-Scaling Test Results 30

5 Weak-Scaling Test Analysis Based on Gustafson’s law: Speedup = f + (1 − f ) × N Serial proportion = 86%, Parallelizable proportion = 14% 31

5 Strong-Scaling Test Results Time required to solve each PDE in each time step 32

5 Strong-Scaling Test Analysis Based on Amdahl’s law: 1 Speedup = f + 1 − f N Serial proportion = 52%, Parallelizable proportion = 48% 33

5 Conclusion ◮ A quantitative mathematical model and its corresponding computational model to assess the degradation behavior of biodegradable materials ◮ Using level set method to track the moving corrosion front during degradation ◮ Once fully validated, the model will be an important tool to find the right design and properties of the metallic biomaterials and implants 34

Thank you for your attention mojtaba.barzegari@kuleuven.be