AN INTEGRATED APPROACH TO MODELING AND MITIGATING SOFC FAILURE - - PowerPoint PPT Presentation

AN INTEGRATED APPROACH TO MODELING AND MITIGATING SOFC FAILURE

Andrei Fedorov, Comas Haynes, Jianmin Qu Georgia Institute of Technology DE-AC26-02NT41571 Program Managers: Travis Shultz National Energy Technology Laboratory

Outline

• First Order Failure Criteria for SOFC PEN Structure
• Creep Modeling of YSZ/Ni Cermet
• Fracture Mechanics Analysis Tool
• Thermal Transient Modeling

First Order Failure Criteria for SOFC PEN Structure

• Objectives
• Local Failure Criteria

– Failure Modes – Strength Failure Criteria – Fracture Failure Criteria

• Global Failure Criteria
• Analyses for Various Crack Cases
• Conclusion

Objectives

Develop first-order failure criteria to be used for the initial design, material selection and

• ptimization against thermomechanical failure
• f the PEN structure in high temperature

SOFCs.

Failure Modes

Material Characteristics

• Static Strength
• Fracture Toughness
• Fatigue Strength

Does a material contain flaws above certain threshold value? No -> Failure is strength- controlled Yes -> Failure is fracture toughness-controlled

Strength-Based Failure Theory

Failure occurs when

1 2 3

( , , )

f

f σ σ σ σ σ = =

where Effective Stress

1 2 3

( , , ) f σ σ σ σ =

1 2 3

, , σ σ σ

Principle Stresses

f

σ

Material Strength

Fracture-Based Failure Theory

Fracture occurs when

2 2 2 2

1 1

III I II c

K G K K G E ν ν ⎛ ⎞ − = + + = ⎜ ⎟ − ⎝ ⎠

where Energy Release Rate G Stress Intensity Factors

I

K

II

K

III

K

Gc Fracture Toughness

YSZ Electrolyte

Maximum Normal Stress Criterion

f

σ σ =

{ }

1 2 3 1 2 3

( , , ) max , , f σ σ σ σ σ σ σ = =

100 ~ 300 MPa

f

σ =

2 2 2 2

1 1

III I II c

K G K K G E ν ν ⎛ ⎞ − = + + = ⎜ ⎟ − ⎝ ⎠

Fracture Criterion

2

7.8 13.7 J m

c

G = ฀

YSZ/Ni Cermet

Von Mises Criterion (elevated temp)

f

σ σ =

( ) ( )

( )

2 2 2 2 2 2

2 2 2

x y y z z x xy yz zx

σ σ σ σ σ σ σ τ τ τ = − + − + − + + +

Maximum Normal Stress Criterion

f

σ σ =

{ }

1 2 3 1 2 3

( , , ) max , , f σ σ σ σ σ σ σ = =

( )

1 (1 )

Ni f YSZ YSZ YSZ Void YSZ Ni

E V V V E σ σ ν ⎡ ⎤ = + − − ⎢ ⎥ − ⎣ ⎦

100 ~ 300MPa

YSZ

σ = = YSZ tensile strength = Ni Young's modulus

Ni

E

= YSZ Young's modulus

YSZ

E

Ni

ν = Ni Poisson's ratio

YSZ

V = YSZ Volume fraction

Void

V

= Void Volume fraction

Global Failure Criteria

Processing Warpage Stack Assembly Fracture

R

1/ R ρ =

c

W W <

c

ρ ρ < W

L

Warpage Criterion Curvature Criterion

Implementation

Based on material/geometry parameters to compute Wc and ρc Measure W or ρ of each cell after sintering Compare the measured W with Wc or ρ with ρc

R

1/ R ρ =

L

W

Crack Types

A – crack in the cathode B – crack in the anode C – delamination crack between the cathode and electrolyte D – delamination crack between the anode and the electrolyte E – blister crack on the anode/electrolyte interface F – crack in the electrolyte

electrolyte cathode anode A F G E electrolyte cathode anode A E C F B

L

ha he hc

D electrolyte cathode anode A F G E electrolyte cathode anode A E C F B

L

ha he hc

D

• Max. Allowable Warpage

Gc = fracture toughness he = electrolyte thickness Ee = modulus of electrolyte

c c e e e

W G L Y L h E h ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠

C rack A

1/ 2 1 / 2 2 2 3 2 2 4 2 3 3

16 (1 ) 2 h E a Y t H E Q α ν

−

⎡ ⎤ ⎛ ⎞ ⎛ ⎞ ∆ ⎛ ⎞ ⎢ ⎥ = + − ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ − ⎝ ⎠ ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦

C rack C

1/ 2 1 / 2 2 3 2 2 3 2 3 2 1 2 2 1 3 3 3

4 ( ) 16 16 c F h F h E Y Q Q h c α

− − −

⎛ ⎞ ⎛ ⎞ ∆ = + ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠

C rack D

1/ 2 1 / 2 2 2 3 1 3 2 2 2

16 (1 ) 2 h h h Y aE Q α π ν

−

⎡ ⎤ ⎛ ⎞ ⎛ ⎞ − ∆ ⎛ ⎞ ⎢ ⎥ = + ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ − ⎝ ⎠ ⎢ ⎥ ⎝ ⎠ ⎝ ⎠ ⎣ ⎦

C rack E

1 / 2 1/ 2 2 3 2 2 2 2 1 2 2 1 3 2

4 ( ) 16 16

ce ce ce ce

c F h F h E Y Q Q h h c α ρ

− − −

⎛ ⎞ ⎛ ⎞ ∆ = + ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠

C rack F

2 2 2 2 1 2 1

1 4

c

Q h h E P P Y P G α ⎛ ⎞ = + ⎜ ⎟ ∆ ⎝ ⎠

Implementation

Definition of Variables: See our monthly report or e-mail jianmin.qu@me.gatech.ed u Basic Assumptions: Linear elastic fracture mechanics Implementation: A FORTRAN code Material Properties Needed: Elastic moduli Coefficient of thermal expansion Fracture toughness Other Parameters needed: Layer thickness Warpage (curvature)

Materials Properties

Young’s Modulus (GPa) Poisson’s Ratio CTE(10-6/oC) Thickness Cathode 90 0.3 11.7 75 Electrolyte 200 0.3 10.8 15 Anode 96 0.3 11.2 500

( ) m µ

Considering sintering process, the set of materials in table will result in

– tensile stress in cathode; – compressive stress in electrolyte; – compressive stress in anode;

Average Stress in Cathode

• 40
• 30
• 20
• 10

10 20 30 40 50 60 20 40 60 80 100 Time Stress (MPa)

Cathode 800oC 800oC 20oC 800oC 20oC cooling heating NiO reduction cooling

Average Stress in Electrolyte

• 180
• 140
• 100
• 60
• 20

20 20 40 60 80 100 Time Stress (MPa)

Electrolyte

800oC 800oC 20oC 800oC 20oC cooling heating NiO reduction cooling Average Stress in Anode

• 6
• 4
• 2

2 4 6 8 20 40 60 80 100 Time Stress (MPa)

Anode

800oC 800oC 20oC 800oC 20oC cooling heating NiO reduction cooling

electrolyte cathode anode A F G E electrolyte cathode anode A E C F B

L

ha he hc

D

• Stress free at 800oC
• No-creep
• NiO reduction results in 0.1%
• vol. shrinkage

Average in-Plane Stress in the PEN Layers

Numerical Examples of Max. Allowable Warpage

Gc (J/m2)

5 10 15 20 W (mm) 2 4 6 10 Crack C Crack D Crack A a/h = 0.01 a/h3 = 0.05 a/h3 = 0.1 8

3

electrolyte cathode anode A F G E electrolyte cathode anode A E C F B

L

h1 h2 h3

D

L = 10 cm

Crack C is the limiting factor, unless crack A is larger than 5% of the cathode thickness.

Statistical Consideration

Failure theories have the following form:

Failure occurs when Σ > Σf

where Σ is the “stress”(e.g., max. normal stress, Mises stress,

• r SIFs, max. warpage, etc. ) and Σf is the“strength”(e.g., yield

strength, fracture toughness, etc.) Both Σ and Σf can be random variables with certain distributions, such normal distribution, Weibull distribution, etc.

Assume: g(σ) = distribution of stress; gf(σ) = distribution of strength The probability of failure at a given stress σ is

( )

f

g x dx

σ −∞

∫

The probability of failure for a given stress distribution g(σ) is

( ) ( )

f f

p g g x dx d

σ

σ σ

∞ −∞ −∞

⎡ ⎤ = ⎢ ⎥ ⎣ ⎦

∫ ∫

Example (Normal Distributions)

2

1 1 ( ) exp 2 2

f f f f

g s s σ σ σ π ⎡ ⎤ ⎛ ⎞ − ⎢ ⎥ = − ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦

Strength distribution

2

1 1 ( ) exp 2 2 g s s σ σ σ π ⎡ ⎤ − ⎛ ⎞ = − ⎢ ⎥ ⎜ ⎟ ⎝ ⎠ ⎢ ⎥ ⎣ ⎦

Stress distribution s = Standard deviation = Mean value

σ

( ) 1 g d σ σ

∞ −∞

=

∫

f

σ σ

2

1 Exp Erfc 2 2 2 2

f f f

p d s s s σ σ σ σ σ π

∞ −∞

⎡ ⎤ ⎡ ⎤ − − ⎛ ⎞ = − ⎢ ⎥ ⎢ ⎥ ⎜ ⎟ ⎝ ⎠ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦

∫

Failure Probability

f

σ σ Factor of Safety

f f f

s s σ σ = Deviation

f

p Failure Probability 1.0 any value 0.5 2.0 0.2

2

3.8 10− × 5.0 0.2

4

2.3 10− × 10.0 0.2

5

7.3 10− × 1.5 0.1

3

9.2 10− × 2.0 0.1

4

2.0 10− × 3.0 0.1

6

1.2 10− × 4.0 0.1

8

5.7 10− × 1.5 0.05

6

1.2 10− × 2 0.05

13

7.7 10− × 1.5 0.02

32

2.3 10− ×

1 2 3 4 0.1 0.2 0.3 0.4 0.5 0.025 0.05 0.075 0.1 1 2 3

f

σ σ

f

s σ

f

p

f

σ σ

f

s s =

Summary

• First Order Failure Criteria
• Local and global failure criterion were established.

These criterion may be easily used to aid the initial design, material selection and optimization of SOFCs.

• Using the local failure criteria, the user can predict

(estimate) the potential material failure

• Using the global failure criteria, the user can predict

whether a cell can survive the stacking assembly process

A Numerical Simulation Tool for Fracture Analysis in Solid Oxide Fuel Cells

( )

1 2

applied stress

i I II

K iK FL

ε −

+ ≡ × K =

KI

• pening

KII shear KIII

• ut of plane

Material 2 Material 1 Crack

Significance of SIFs

1. Will the crack grow? 2. In what direction? (What is mode mixity?)

( )

2 * 2 *

1 cosh 2

III ic

K G E πε µ + = KK

1 Im[

] tan Re[ ]

i i

L L

ε ε

ψ

− ⎡

⎤ = ⎢ ⎥ ⎣ ⎦ K K

Temperature

Computing Fracture Parameters Using Volume Integrals

int , , , , , kj aux aux aux aux ij ij k ij j ik ij i j k ii k j

P u u x σ ε σ σ ασ θ = − − − ∂

int aux aux aux i i jk mn mn jk ik ik j j

u u P x x σ ε δ σ σ ∂ ∂ = − − ∂ ∂

Curvilinear Strain Energy Stress and spatial derivatives of displacements

int int int j kj jk j V k k

q P I G P q dV x x ⎛ ⎞ ∂ ∂ = = − + ⎜ ⎟ ⎜ ⎟ ∂ ∂ ⎝ ⎠

∫

Virtual crack growth (Q)

dV

x1 y1

Crack u, σ, and ε from FEM uaux, σaux, and εaux analytical dV

x1 y1

Crack u, σ, and ε from FEM uaux, σaux, and εaux analytical

X1 X3 z1 x1 a Crack Edge Volume Boundary

Q

X1 X3 z1 x1 a Crack Edge Volume Boundary

Q

• Pointwise Value
• To find KI by setting
• KI

aux=1

• KII

aux=KIII aux=0

* 2 *

2 1 ( ) cosh ( )

aux aux aux I I II II III III

I s K K K K K K E πε µ ⎡ ⎤ = + + ⎣ ⎦

( )

* 2

( ) cosh 2

I

I s K E πε =

( ) ( )

c

L

I I s a s ds = ∆

∫

A Pe nny-Shape d Cr ac k on E le c tr

• lyte / Anode Inte r

fac e

Half Model - Interface Penny Crack

w=15a a h=30a e= a/14 18 layers e = a/14 18 layers Anode Electrolyte

Mesh footprint near crack

θ Crack Tip z y x

Half Model - Interface Penny Crack

w=15a a h=30a e= a/14 18 layers e = a/14 18 layers Anode Electrolyte

Mesh footprint near crack

θ Crack Tip z y x

Temperature gradient parallel to crack plane

T e mpe r atur e Gr adie nt Par alle l to the Cr ac k Plane

Tmax T=0ºC

ANSYS Model FMA Volume

x

Electrolyte Anode

σy (MPa)

K

I, K II, K III Var

iation Along Cr ac k F r

• nt
• 0.015
• 0.01
• 0.005

0.005 0.01 0.015 30 60 90 120 150 180 Angle Normalized Value KI KII KIII

σxy (MPa)

Summary

• - Fracture Mechanics Analysis Tool
• Based on volume integral (requires less mesh density)
• Written in MatLab language (run on both Window and Unix)
• Add-on to any commercial FEM codes (requires less

processing time) Fracture Mechanics Analysis Tool: Capabilities:

• Calculate energy release rate and individual stress intensity

factors

• 2D and 3D planar cracks of arbitrary shapes
• Homogeneous and interfacial cracks
• Arbitrary mechanical and thermal loading

Transient Heat Transfer Analysis: Convective-Conductive Heating of SOFC

SOFC unit cell Transient Thermal Modeling

Key Question/Focus: Provide model-based design tool(s) to assess how quickly a cell/stack cane be heated without excessive (damaging) thermomechanical gradients?

Design/Model Outputs:

• total time required for heating
• max temperature spatial gradient
• max temperature time-derivative

Model Inputs:

• size of components; thickness of

layers

• thermal properties of components
• boundary conditions; heating

strategies Multi-Level Methodology:

• 3-D CFD modeling (e.g. FLUENT)
• Reduced order numerical modeling
• Simplified order analytical modeling

Current Collector Current Collector Electrolyte Cathode Anode Air Channel Fuel Channel

y x z

Current Collector Current Collector Electrolyte Cathode Anode Air Channel Fuel Channel Current Collector Current Collector Electrolyte Cathode Anode Air Channel Fuel Channel

y x z y x z

Simplified Analytical Model/Design Tool: Key Ideas & Assumptions

Heating by hot air supplied at prescribed time-dependent inlet temperature

• 1-D temperature profile in each component Tlayer = f (z,t) only
• Constant velocity

plug flow in channel

• Constant properties
• Radiation neglected

(to be included later)

• Adiabatic boundaries

(no heat losses)

L Hot Air @ const V Fuel Channel Interconnect z Anode Electrolyte Cathode Interconnect Insulated Boundaries

1st order 1-D Purely Convective Heating Model

Key Assumptions:

• thermally thin cell materials

(i.e. no energy storage)

• thermal equilibrium between

air and channel walls

z c.v u , T hot T w all = T air z c.v u , T hot T w all = T air

B.C. & I.C.

( ) ( )

0, ( , 0)

• T z

t f t T z t T = = = =

T T u t z δ δ δ δ + =

Governing Equation ( )

for , for

• T

z ut T z t z f t z ut u > ⎧ ⎪ = ⎨ ⎛ ⎞ − ≤ ⎜ ⎟ ⎪ ⎝ ⎠ ⎩

Closed-form analytical solution:

2nd Order Convective-Conductive Heating Model

( )

( )

( ) ( ) ( )

( )

( ) ( ) ( )

( )

( )

1 1 1 1 1 1 1 1 1 1 1

1 2 2 2 2 2 2

Air InterConne Channe ct : Cat l (gas) hode: :

IC C p g IC IC IC IC C p g C g IC g g IC IC IC IC C IC C g C g g g g g g IC p C g C C C C

T T T T T P c A kA hP t z R c A u kA hP hP t z z P c A k T T T A hP t T T T T T z T T T T ρ ρ ρ

− − − − − −

∂ ∂ = + − − − ∂ ∂ ∂ ∂ ∂ ⎡ ⎤ + = − − − − ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦ ∂ ∂ = + − + ∂ ∂

( )

( )

( )

( ) ( ) ( )

( )

( ) ( )

( ) ( ) ( )

( )

1 1 1 2 2 2

2 2 2 2 2

Electrolyte: Anode: FuelChannel:

C C E IC C C E C E E E E A p C E E IC C C A E E C E E A A IC A A E A p E A f A f A A IC A A E A A IC fuel p fuel fue E l

P R R P T T P c A kA T T T T t z R R P T T P c A kA T T hP T T T T t z T R T R T c A kA T t T ρ ρ ρ

− − − − − − − − − − − − −

− − − ∂ ∂ = + − − − ∂ ∂ ∂ ∂ = + − − − − − ∂ ∂ ∂ ∂ = ∂

( ) ( )

2 2

2 fuel f A f A f IC f IC

T hP T T hP T T z

− −

+ − + − ∂

Applying thermal equilibrium between flow channels and components; model reduces to a single equation dependent only on effective Peclet number and inlet temperature function!

2 2

1 T T T t z z Pe ∂ ∂ ∂ + = ∂ ∂ ∂ ( ) ( ) ( ) ( ) ( )

. .& . .: 0, 0, ; 1, 0; ,0 1 T T B C I C T t t t Pe t z z z F T ∂ ∂ − = = ∂ ∂ advection conductio

• f thermal energy

n

eff eff

u L Pe α = ≡ = Closed-form analytical solution has been obtained!!!

Results: Comparison of 1st and 2nd order models

Key advantages demonstrated:

• Computationally efficient, analytical models capture key physics of heating process!
• 1st order model is the limiting case of 2nd order model ( Pe large ) provided it is

properly re-scaled. The guidelines for re-scaling have been developed!

1.0E-03 1.0E-02 1.0E-01 1.0E+00 1.0E+01 1.0E+01 1.0E+02 1.0E+0

Thermally thin limit

Pe = 0.1 Pe = 0.5 Pe = 1

Heating Time T z ∂ ∂

K = 0.1 K = 0.5 K = 1 K = 2 K = 0.2 K = 0.05

→ ∞ Pe

10 1 0.1 0.01 10 100 1000 0.01 0.10 1.00 10.00 1 10 100 1000

Pe = 0.1 Pe = 10 Pe = 1

Heating Time T t ∂ ∂

K = 0.1 K = 0.5 K = 1 K = 10 K = 0.05

→∞ Pe

0.01 0.10 1.00 10.00 1 10 100 1000 Pe = 0.1 Pe = 10 Pe = 1

Heating Time T t ∂ ∂

K = 0.1 K = 0.5 K = 1 K = 10 K = 0.05

→∞ Pe

Design Maps: Dimensionless plots of temperature gradient and time-derivative vs. total heating time for various rates of inlet temperature rise (K) and Peclet numbers (Pe).

Summary

• - Convective-Conductive Heating of SOFC

Developed reduced order solutions for transient thermal analysis. Obtained closed-form analytical solutions that provide a relationship between heating rate and the spatial temperature gradient Obtained closed-form analytical solutions that provide a relationship between rate and the temporal temperature gradient

Future Work

Refine and validate the first order failure criteria Develop and implement the global-local computational algorithm in MARC Validate and implement a suitable FEA tool for analysis of fracture failure in the context of various pre-existing flaws within SOFC cells under various operating conditions. Validate and implement a computationally- efficient transient thermal model.