Species Transport

Overview

The species transport module (solver_species/mod_species.f90) enables simulation of dissimilar metal mixing in additive manufacturing. It solves a convection-diffusion equation for concentration and provides composition-dependent material properties through two-way coupling.

Secondary material properties and transport numerics live in solver_species/inputfile/input_param_species.txt and are read by read_species_params() at startup when species_flag=1.

Concentration Convention

Value	Meaning	Properties Source
C = 1	Pure primary (base) material	&material_properties in input_param.txt
C = 0	Pure secondary material	&species_secondary_material / &species_secondary_powder in input_param_species.txt
0 < C < 1	Mixed region	mix(prop1, prop2, C) = prop1*C + prop2*(1-C)

Enabling Species Transport

In input_param.txt:

&output_control ... species_flag=1 /

When species_flag=0 (default), the species module is completely inert — no arrays allocated, no computation, no overhead.

Physics

Transport Equation

\[\frac{\partial (\rho C)}{\partial t} + \nabla \cdot (\rho \vec{u} C) = \nabla \cdot (\Gamma \nabla C)\]

where \(\Gamma = \rho D_m\) is the mass diffusion coefficient, \(D_m = 5 \times 10^{-9}\) m\(^2\)/s.

In solid regions (\(T < T_s\)), \(\Gamma\) is set to a floor value (\(10^{-30}\)) to effectively freeze concentration.

Two-Way Coupling

When species_flag=1, material properties depend on local concentration:

• Thermal: tsolid, tliquid, acpa, acpb, acpl, thconsa, thconsb, thconl
• Mechanical: dens, denl, viscos
• Powder: pden, pcpa, pcpb, pthcona, pthconb

All computed inline via mix() — no extra arrays needed.

Properties that remain scalar (same for both materials):

• beta (thermal expansion), emiss (emissivity)
• hlatnt (latent heat), dgdt (thermal Marangoni)

Solutal Marangoni

Surface tension depends on concentration: \(\gamma = \gamma_0 + \frac{d\gamma}{dT}(T - T_{ref}) + \frac{d\gamma}{dC}(C - C_{ref})\)

The solutal Marangoni stress drives flow from low-C to high-C regions (when dgdc < 0):

\[\tau_C = \frac{d\gamma}{dC} \frac{\partial C}{\partial x}\]

Currently dgdc_const = -0.3 N/m (typical for liquid metal binary alloys).

Initial Conditions

• Substrate (below powder layer): C = 1 (primary material everywhere)
• Powder layer, y < y_mid: C = 1 (primary)
• Powder layer, y ≥ y_mid: C = 0 (secondary)

where y_mid = dimy / 2. This represents a half-and-half powder bed configuration for dissimilar metal testing — substrate is uniform primary, the powder layer is split along the Y-axis.

The IC is built directly on the species-private uniform grid in species_init_concentration_sp(), then projected onto the thermal grid via species_interp_species_to_thermal() so the first thermal step sees a populated concentration field. After init, rebuild_enthalpy_from_temp() (in mod_entot) re-derives every cell's enthalpy from tempPreheat under the local material mix — without this step, secondary-material cells round-trip the primary-mix preheat enthalpy through enthalpy_to_temp and land at non-physical temperatures (e.g., a C=0 cell with primary preheat enthalpy comes back as 247.8 K instead of 293 K under acpa2=0.3441, acpb2=400).

Secondary Material Properties

Read from solver_species/inputfile/input_param_species.txt via namelists &species_secondary_material, &species_secondary_powder, and &species_transport. Defaults below match the historical hardcoded constants:

Property	Primary	Secondary	Unit
dens/dens2	8440	8880	kg/m^3
denl/denl2	7640	7800	kg/m^3
viscos/viscos2	0.007	0.003	Pa*s
tsolid/tsolid2	1563	1728	K
tliquid/tliquid2	1623	1803	K
acpa/acpa2	0.2441	0.3441	J/(kg*K^2)
acpb/acpb2	338.59	400	J/(kg*K)
acpl/acpl2	709.25	800	J/(kg*K)
thconl/thconl2	30.078	120	W/(m*K)
D_m	—	5e-9	m^2/s
dgdc_const	—	-0.3	N/m

Solver Details

• Grid: species runs on its OWN uniform X/Y grid (not the thermal AMR grid). Cell size is amr_dx_fine when adaptive_flag=1, or dimx/ncvx_orig(1) when AMR is off. Z is shared 1:1 with the thermal solver. Built once in species_init_grid().
• Discretization: Power-law convection-diffusion scheme on the species grid
• Solver: Line-by-line TDMA + block correction (solve_species)
• Domain: full species grid (no melt-pool clipping). The fine, decoupled grid is cheap enough that bbox-restriction is unnecessary — correctness first.
• Timing: once per thermal timestep, after the thermal iteration loop exits
• Under-relaxation: urfspecies = 0.7
• Clipping: concentration forced to [0, 1] after each solve
• Transient: uses delt (shared with thermal)

Two-Way Coupling Step Sequence

Each thermal timestep, after the iteration loop:

1. species_compute_thermal_drivers — bilinear-interpolate temp, uVel, vVel, wVel, den, vis from the (variable, AMR-shrunk) thermal grid to the (fixed, uniform) species grid.
2. species_bc — Neumann zero-flux BCs on all 6 faces of the species grid.
3. solve_species — one transport sweep on the species grid (TDMA + block correction).
4. species_interp_species_to_thermal — area-weighted average of concentration_sp back onto the thermal grid so the next step's properties() / source_enthalpy() / enthalpy_to_temp() calls read a fresh, AMR-aware concentration(i,j,k) for their mix() calls.

The thermal-grid concentration array is therefore a derived, frequently-refreshed view of concentration_sp — never modified directly by anyone else.

Performance

Typical overhead with species enabled:

Metric	Impact
Species transport solve	~3-5% of total wall (full-grid TDMA is cheap)
Two-way coupling helpers	~5% (bilinear interp + area-weighted reproject)
mix() overhead in properties()/source/enthalpy_to_temp	~10%
Memory	+60-100 MB depending on amr_dx_fine and domain size

Output

When species_flag=1, the species solver writes its OWN file series to result/<case>/species_results/, decoupled from the thermal solver's output:

File	Format	Description
<case>_species{N}.vts	XML StructuredGrid + b64-zlib	Per-frame snapshot on the species uniform grid
<case>_species.pvd	XML Collection	Time-aware index — open this in ParaView

Per-frame fields:

• Velocity (vector, m/s) — interpolated from the thermal grid
• T (K) — interpolated from the thermal grid
• concentration (-) — mass fraction of primary material [0, 1]

den, vis, and tsolid_field are intentionally NOT in the species output: den/vis are already in the thermal vtkmov*.vts, and tsolid_field = mix(tsolid, tsolid2, concentration) is a pure ParaView Calculator expression away (tsolid * concentration + tsolid2 * (1 - concentration)).