Complete implementation of the finite element model for gravity-induced grain settling and shape distortion during liquid phase sintering of tungsten heavy alloys.
Based on: Alvarado-Contreras et al. / Acta Materialia 65 (2014) 176-184
"A continuum approach for modeling gravitational effects on grain settling and shape distortion during liquid phase sintering of tungsten heavy alloys"
This code implements ALL equations and parameters from the paper without any simplifications:
-
material_parameters.py: All material parameters from Table 1
- Average grain radius: 10 μm
- Surface tension: 2.5 J m⁻²
- Gravity: 9.81 m s⁻²
- Friction coefficient: 0.80
- Densities, viscosities, critical volume fractions
- Arrhenius parameters from Figure 2
-
constitutive_equations.py: Constitutive formulation
- Equation (2): Stress-strain rate relationship
- Equation (3): Strain rate tensor
- Equation (4): Viscosity definitions
- Equation (5): Normalized viscosities
- Equation (6): Effective sintering stress
- Equation (7): Viscosity estimation from dilatometry
- Equation (8): Arrhenius viscosity-temperature relationship
- Equation (9): Activation energy
- Equation (10): Pre-exponential term with volume fraction dependence
- Equation (19): Shear rate magnitude
-
mass_balance.py: Mass conservation equations
- Equation (12): Mass balance for solid and liquid phases (fully implemented)
- Equation (13): Densification rate (COMPLETE FORM - all parameters required, no simplifications)
- Equation (14): Porosity kinetics (complete implementation uses Equation 13)
- Equation (15): Local expansion rate due to grain flux
- Equation (16): Velocity decomposition
- Equation (17): Solid volume fraction equation (Kynch-like)
-
grain_migration.py: Grain migration and flux
- Equation (18): Concentration and viscosity gradient fluxes
- Equation (19): Shear rate magnitude
- Equation (20): Stokes' settling velocity
- Equation (21): Settling flux
- Equation (22): Hindered function
- Equation (23): Total flux density
-
finite_element_model.py: Finite element implementation
- Axisymmetric domain
- Second-order isoparametric quadrilateral elements
- Equation (24): Boundary conditions (friction, fixed displacement)
- Equation (25): Surface tension
- Zero grain flux boundary condition at specimen surface (from paper)
-
temperature_profile.py: Heating schedule
- Temperature range: 900-1500 °C
- Heating rate: 15 °C/min
- Holding time: 30 min at 1500 °C
-
main_simulation.py: Complete simulation program
- Integrates all modules
- Time-stepping algorithm
- Results storage and visualization
σ_ij = 2η(ε_ij - (1/3)εδ_ij) + ζεδ_ij + P_Lδ_ij
η = φ(θ)η_0(φ, T)
ζ = 2ψ(θ)η_0(φ, T)
φ = (1 - θ)²
ψ = (2/3)(1 - θ)³/θ
P_L = 3α(1 - θ)²/r_0
ln η_0 = ln A + B/(RT) + C/(RT)²
A = A_88W(1 - φ/φ_c)^(-mφ_c)
∂(1-φ)/∂t + ∇·((1-φ)u) = -∇·J
∂θ/∂t = -∇·[(1-θ)u] + (Δρ/(ρ_S^0 ρ_L^0)) * ∇·[J * (φ - θ)] (COMPLETE FORM)
∂φ/∂t + ∇·(φu) = -∇·J
J_i^c = -ρ_s K_c a² ∇_i(γ̇φ)
J_i^v = -ρ_s K_η γ̇φ² (a²/η_0) ∇_i(η_0)
J_i^s = ρ_s f_H u̇_i^∞ φ
J_i = J_i^c + J_i^v + J_i^s
u_z = 0, σ_r = μσ_z (bottom surface)
S_i^n = -ακn_i (free surface)
- Average grain radius, a: 10 μm
- Surface tension, σ: 2.5 J m⁻²
- Gravity acceleration, g: 9.81 m s⁻¹
- Coulomb frictional factor, μ: 0.80
- Tungsten density, ρ_s: 19.30 g cm⁻³
- Liquid phase density, ρ_L: 8.57 g cm⁻³
- Initial solid volume fraction, Φ: 0.764
- Initial porosity, θ₀: 0.34
- Grain size, r₀: 10 μm
- Liquid phase viscosity, η_L: 0.25 Pa s
- Critical volume fraction, Φ_c: 0.80
- Viscosity exponent, m: 2.5
- Hindered function exponent, n: 3.25
- Ln A = 19.195
- B/R = 21334.0 K
- C/R² = 0 (solid-state)
- C/R² = -2.5×10⁵ (liquid-phase)
- Transition temperature, T_t = 1460 °C
from main_simulation import SinteringSimulation
# Create simulation for 88W alloy
sim = SinteringSimulation(W_content=88, aspect_ratio=1.0, D0=10e-3)
# Run simulation
results = sim.run_simulation(dt=0.1, save_interval=1.0)
# Plot results
sim.plot_results()- Sample geometry: Cylindrical, 10 mm diameter × 10 mm height (or variable aspect ratio)
- Temperature: 900-1500 °C
- Heating rate: 15 °C/min
- Holding: 30 min at 1500 °C
- Tungsten content: 83, 88, or 93 wt%
The simulation produces:
- Porosity evolution vs. temperature
- Volume fraction evolution vs. temperature
- Final geometry (distorted shape)
- Top-bottom volume fraction gradients
- All state variables at each time step
- numpy
- matplotlib
- scipy (for advanced FE operations, if needed)
This implementation follows the paper exactly with no simplifications. All equations, parameters, and boundary conditions are implemented as specified in the original publication.
- ✅ All 25 equations from the paper are implemented
- ✅ All parameters from Table 1 are exactly as specified
- ✅ Zero grain flux boundary condition at specimen surface (as mentioned in paper)
- ✅ Proper axisymmetric formulation including u_r/r terms in divergence
- ✅ Complete mass balance for solid and liquid phases (Equation 12)
- ✅ Velocity decomposition (Equation 16) for sintering and grain transport
- ✅ Viscosity estimation from dilatometry (Equation 7) for calibration
- ✅ Momentum conservation equation properly formulated (div(σ) + f = 0)
- The code uses proper axisymmetric divergence:
div = d(v_r)/dr + v_r/r + d(v_z)/dz - All boundary conditions from the paper are implemented, including the zero grain flux condition
- The Arrhenius viscosity model includes both solid-state and liquid-phase regions
- Grain migration includes all three flux components: concentration, viscosity gradient, and settling
This project is licensed under the MIT License - see the LICENSE file for details.
Contributions are welcome! Please feel free to submit issues or pull requests.