diff --git a/src/ProcessSimulator.jl b/src/ProcessSimulator.jl
index 2ae7646..d84ae48 100644
--- a/src/ProcessSimulator.jl
+++ b/src/ProcessSimulator.jl
@@ -14,6 +14,10 @@ include("fluid_handling/compressors.jl")
 include("fluid_handling/heat_exchangers.jl")
 
 # Reactors
+include("reactors/types.jl")
+include("reactors/reaction_sets.jl")
 include("reactors/CSTR.jl")
+include("reactors/PFR.jl")
+include("reactors/GibbsReactor.jl")
 
 end
diff --git a/src/reactors/GibbsReactor.jl b/src/reactors/GibbsReactor.jl
new file mode 100644
index 0000000..965af85
--- /dev/null
+++ b/src/reactors/GibbsReactor.jl
@@ -0,0 +1,209 @@
+"""
+    GibbsReactor - Equilibrium reactor based on Gibbs free energy minimization
+
+The GibbsReactor finds the equilibrium composition by minimizing the total
+Gibbs free energy subject to elemental balance constraints. This approach
+does not require specification of reaction stoichiometry - equilibrium is
+determined purely from thermodynamic data.
+
+Returns an `OptimizationSystem` for solving the equilibrium problem.
+"""
+
+"""
+    GibbsReactor(ms::MaterialSource; name, element_matrix, feed_elements)
+
+Create a Gibbs equilibrium reactor that minimizes Gibbs free energy.
+
+# Arguments
+- `ms::MaterialSource`: Material source with thermodynamic property functions
+- `name`: System name (required)
+- `element_matrix::Matrix{Float64}`: Element composition matrix (N_elements × N_components)
+  where entry [i,j] is the number of atoms of element i in component j
+- `feed_elements::Vector{Float64}`: Total moles of each element in the feed
+
+# Returns
+An `OptimizationSystem` that, when solved, gives the equilibrium composition.
+
+# Notes
+The Gibbs free energy is computed as G = H - T*S where H and S are obtained
+from the MaterialSource enthalpy and entropy functions. The optimization
+problem minimizes total G subject to:
+1. Element balance: A * n = b (conservation of atoms)
+2. Non-negativity: n ≥ 0 (no negative mole amounts)
+
+# Example
+```julia
+# Define element matrix for H2O, H2, O2 system
+# Rows: [H, O], Columns: [H2O, H2, O2]
+elem_matrix = [2.0 2.0 0.0;   # H atoms
+               1.0 0.0 2.0]   # O atoms
+
+# Feed composition (in terms of elements)
+feed_elem = [4.0, 2.0]  # 4 mol H, 2 mol O
+
+@named gibbs = GibbsReactor(ms; element_matrix=elem_matrix, feed_elements=feed_elem)
+```
+"""
+@component function GibbsReactor(ms::MaterialSource;
+        name,
+        element_matrix::Matrix{Float64},
+        feed_elements::Vector{Float64},
+        T_specified::Union{Nothing, Float64} = nothing,
+        p_specified::Float64 = 101325.0)
+    N_c = ms.N_c
+    N_elem = size(element_matrix, 1)
+
+    # Validate dimensions
+    if size(element_matrix, 2) != N_c
+        throw(ArgumentError(
+            "Element matrix must have $N_c columns (number of components), " *
+            "got $(size(element_matrix, 2))"))
+    end
+    if length(feed_elements) != N_elem
+        throw(ArgumentError(
+            "feed_elements must have $N_elem entries (number of elements), " *
+            "got $(length(feed_elements))"))
+    end
+
+    # Parameters
+    pars = @parameters begin
+        p = p_specified, [description = "System pressure (Pa)"]
+        T_spec = isnothing(T_specified) ? 298.15 : T_specified,
+        [description = "Specified temperature (K), used if isothermal"]
+        A[1:N_elem, 1:N_c] = element_matrix, [description = "Element composition matrix"]
+        b[1:N_elem] = feed_elements, [description = "Feed element amounts (mol)"]
+    end
+
+    # Variables for optimization
+    vars = @variables begin
+        (n(t))[1:N_c], [description = "Molar amounts of each component (mol)", bounds = (0, Inf)]
+        n_total(t), [description = "Total moles", bounds = (0, Inf)]
+        (x(t))[1:N_c], [description = "Mole fractions", bounds = (0, 1)]
+        T(t), [description = "Temperature (K)", bounds = (200, 5000)]
+        G_total(t), [description = "Total Gibbs free energy (J)"]
+        H_total(t), [description = "Total enthalpy (J)"]
+        S_total(t), [description = "Total entropy (J/K)"]
+    end
+
+    # Compute molar density for property calculations
+    ρ_func = (p_val, T_val, x_val) -> ms.molar_density(p_val, T_val, x_val; phase = "unknown")
+
+    # Build equations
+    eqs = Equation[]
+
+    # Total moles
+    push!(eqs, n_total ~ sum(collect(n)))
+
+    # Mole fractions (with protection against division by zero)
+    for i in 1:N_c
+        push!(eqs, x[i] ~ n[i] / n_total)
+    end
+
+    # Element balance constraints: A * n = b
+    for i in 1:N_elem
+        push!(eqs, sum(scalarize(A[i, :] .* n)) ~ b[i])
+    end
+
+    # Thermodynamic properties
+    # Enthalpy: H = sum(n_i * h_i) where h_i is molar enthalpy
+    push!(eqs, H_total ~ n_total * ms.VT_enthalpy(ρ_func(p, T, collect(x)), T, collect(x)))
+
+    # Entropy: S = sum(n_i * s_i) where s_i is molar entropy
+    push!(eqs, S_total ~ n_total * ms.VT_entropy(ρ_func(p, T, collect(x)), T, collect(x)))
+
+    # Gibbs free energy: G = H - T*S
+    push!(eqs, G_total ~ H_total - T * S_total)
+
+    # Temperature specification (if isothermal)
+    if !isnothing(T_specified)
+        push!(eqs, T ~ T_spec)
+    end
+
+    # Create OptimizationSystem
+    # Objective: minimize G_total
+    # Subject to: element balance constraints (included in eqs)
+    return OptimizationSystem(
+        G_total,  # Objective function to minimize
+        eqs,
+        t,
+        collect(Iterators.flatten(vars)),
+        collect(Iterators.flatten(pars));
+        name
+    )
+end
+
+"""
+    SteadyStateGibbsReactor(ms::MaterialSource; name, kwargs...)
+
+Create a steady-state Gibbs reactor with inlet and outlet connectors.
+This wraps the GibbsReactor optimization in a component that can be
+connected to other process units.
+
+# Arguments
+- `ms::MaterialSource`: Material source with thermodynamic property functions
+- `name`: System name (required)
+- `element_matrix::Matrix{Float64}`: Element composition matrix
+- `feed_elements::Vector{Float64}`: Feed element amounts (will be computed from inlet)
+
+# Returns
+An ODESystem component with material connectors.
+"""
+@component function SteadyStateGibbsReactor(ms::MaterialSource;
+        name,
+        element_matrix::Matrix{Float64})
+    N_c = ms.N_c
+    N_elem = size(element_matrix, 1)
+
+    # Connectors
+    @named In = MaterialConnector(ms)
+    @named Out = MaterialConnector(ms)
+
+    # Parameters
+    pars = @parameters begin
+        A[1:N_elem, 1:N_c] = element_matrix, [description = "Element composition matrix"]
+    end
+
+    # Variables
+    vars = @variables begin
+        (n_out(t))[1:N_c], [description = "Equilibrium molar amounts (mol)"]
+        n_total_out(t), [description = "Total outlet moles"]
+        (elem_feed(t))[1:N_elem], [description = "Element amounts in feed (mol)"]
+        G(t), [description = "Gibbs free energy at equilibrium (J)"]
+    end
+
+    eqs = Equation[]
+
+    # Compute element amounts from inlet composition
+    for i in 1:N_elem
+        push!(eqs, elem_feed[i] ~ sum(scalarize(A[i, :] .* In.xᵢ .* In.n)))
+    end
+
+    # Element balance at outlet (steady state: in = out)
+    for i in 1:N_elem
+        push!(eqs, elem_feed[i] ~ sum(scalarize(A[i, :] .* n_out)))
+    end
+
+    # Total outlet moles
+    push!(eqs, n_total_out ~ sum(collect(n_out)))
+
+    # Outlet connector equations
+    push!(eqs, Out.n ~ -n_total_out)  # Negative for outflow
+    push!(eqs, Out.p ~ In.p)
+    push!(eqs, Out.T ~ In.T)  # Isothermal assumption for this simplified version
+    for i in 1:N_c
+        push!(eqs, Out.xᵢ[i] ~ n_out[i] / n_total_out)
+    end
+
+    # Flow balance
+    push!(eqs, In.n + Out.n ~ 0)
+
+    # Gibbs energy (for output/monitoring)
+    push!(eqs, G ~ ms.VT_enthalpy(Out.ϱ, Out.T, collect(Out.xᵢ)) * n_total_out -
+                   Out.T * ms.VT_entropy(Out.ϱ, Out.T, collect(Out.xᵢ)) * n_total_out)
+
+    return ODESystem(eqs, t, collect(Iterators.flatten(vars)),
+        collect(Iterators.flatten(pars));
+        name, systems = [In, Out])
+end
+
+export GibbsReactor, SteadyStateGibbsReactor
diff --git a/src/reactors/PFR.jl b/src/reactors/PFR.jl
new file mode 100644
index 0000000..27ea0df
--- /dev/null
+++ b/src/reactors/PFR.jl
@@ -0,0 +1,311 @@
+"""
+    PFR - Plug Flow Reactor
+
+The Plug Flow Reactor (PFR) models a tubular reactor with no axial mixing.
+The composition varies along the reactor length as reactions proceed.
+This implementation uses spatial discretization to convert the PDEs into
+a system of ODEs.
+
+Key assumptions:
+- No radial gradients (1D model)
+- No axial dispersion (plug flow)
+- Steady-state or dynamic operation
+- Constant cross-sectional area (can be extended)
+"""
+
+"""
+    PFR(ms::MaterialSource; name, L, A, N_nodes, i_reacts, adiabatic)
+
+Create a Plug Flow Reactor component.
+
+# Arguments
+- `ms::MaterialSource`: Material source with thermodynamic properties and reactions
+- `name`: System name
+- `L::Float64`: Reactor length (m)
+- `A::Float64`: Cross-sectional area (m²)
+- `N_nodes::Int`: Number of discretization nodes along length
+- `i_reacts`: Indices of reactions to include (default: all)
+- `adiabatic::Bool`: If true, no heat transfer; if false, Q must be specified (default: true)
+
+# Returns
+An `ODESystem` representing the discretized PFR.
+
+# Notes
+The PFR is discretized into N_nodes spatial nodes. At each node, mass and
+energy balances are solved. The spatial coordinate z varies from 0 to L.
+
+For steady-state operation, use with `structural_simplify` and set
+time derivatives to zero.
+
+# Example
+```julia
+@named pfr = PFR(ms; L=2.0, A=0.01, N_nodes=20)
+```
+"""
+@component function PFR(ms::MaterialSource;
+        name,
+        L::Float64,
+        A::Float64,
+        N_nodes::Int = 20,
+        i_reacts = 1:length(ms.reaction),
+        adiabatic::Bool = true)
+    # Validate
+    N_nodes >= 2 || throw(ArgumentError("N_nodes must be at least 2"))
+    L > 0 || throw(ArgumentError("Reactor length L must be positive"))
+    A > 0 || throw(ArgumentError("Cross-sectional area A must be positive"))
+
+    N_c = ms.N_c
+    N_r = length(ms.reaction[i_reacts])
+    Δz = L / (N_nodes - 1)  # Node spacing
+    V_node = A * Δz  # Volume per node
+
+    # Connectors
+    @named In = MaterialConnector(ms)
+    @named Out = MaterialConnector(ms)
+
+    # Parameters
+    pars = @parameters begin
+        L_reactor = L, [description = "Reactor length (m)"]
+        A_cross = A, [description = "Cross-sectional area (m²)"]
+        dz = Δz, [description = "Node spacing (m)"]
+        V_seg = V_node, [description = "Volume per segment (m³)"]
+    end
+
+    # Variables at each node
+    vars = @variables begin
+        # State variables at each node
+        (T_node(t))[1:N_nodes], [description = "Temperature at each node (K)"]
+        (p_node(t))[1:N_nodes], [description = "Pressure at each node (Pa)"]
+        (xᵢ_node(t))[1:N_nodes, 1:N_c], [description = "Mole fractions at each node"]
+        (ϱ_node(t))[1:N_nodes], [description = "Molar density at each node (mol/m³)"]
+        (n_node(t))[1:N_nodes], [description = "Molar flow at each node (mol/s)"]
+        (h_node(t))[1:N_nodes], [description = "Molar enthalpy at each node (J/mol)"]
+
+        # Reaction rates at each node
+        (r_node(t))[1:N_nodes, 1:N_r], [description = "Reaction rates at each node (mol/s)"]
+        (R_node(t))[1:N_nodes, 1:N_c], [description = "Species production rates (mol/m³/s)"]
+        (ΔHᵣ_node(t))[1:N_nodes], [description = "Heat of reaction at each node (J/s)"]
+
+        # Outlet variables
+        F_out(t), [description = "Outlet molar flow (mol/s)"]
+        T_out(t), [description = "Outlet temperature (K)"]
+
+        # Heat transfer (if not adiabatic)
+        (Q_node(t))[1:N_nodes], [description = "Heat transfer at each node (J/s)"]
+    end
+
+    eqs = Equation[]
+
+    # Boundary conditions at inlet (node 1)
+    push!(eqs, T_node[1] ~ In.T)
+    push!(eqs, p_node[1] ~ In.p)
+    push!(eqs, n_node[1] ~ In.n)
+    for i in 1:N_c
+        push!(eqs, xᵢ_node[1, i] ~ In.xᵢ[i])
+    end
+
+    # Equations at each node
+    for j in 1:N_nodes
+        # Density from EOS
+        push!(eqs, ϱ_node[j] ~ ms.molar_density(p_node[j], T_node[j],
+            collect(xᵢ_node[j, :]); phase = "unknown"))
+
+        # Enthalpy from EOS
+        push!(eqs, h_node[j] ~ ms.VT_enthalpy(ϱ_node[j], T_node[j], collect(xᵢ_node[j, :])))
+
+        # Reaction rates at this node
+        for (ir, reac) in enumerate(ms.reaction[i_reacts])
+            push!(eqs, r_node[j, ir] ~ reac.r(p_node[j], T_node[j], collect(xᵢ_node[j, :])))
+        end
+
+        # Species production rates
+        for i in 1:N_c
+            push!(eqs, R_node[j, i] ~ sum([r_node[j, ir] * ms.reaction[i_reacts][ir].ν[i]
+                                           for ir in 1:N_r]))
+        end
+
+        # Heat of reaction
+        push!(eqs, ΔHᵣ_node[j] ~ sum([r_node[j, ir] * ms.reaction[i_reacts][ir].Δhᵣ(T_node[j])
+                                      for ir in 1:N_r]) * V_seg)
+
+        # Heat transfer (adiabatic or specified)
+        if adiabatic
+            push!(eqs, Q_node[j] ~ 0.0)
+        end
+    end
+
+    # Spatial derivatives using finite differences (upwind scheme)
+    # dF_i/dz = R_i * A  (component molar flow)
+    # dT/dz depends on energy balance
+    for j in 2:N_nodes
+        # Component balances: d(n*x_i)/dz = R_i * A
+        # Using upwind finite difference: (n*x_i)[j] - (n*x_i)[j-1] = R_i * A * dz
+        for i in 1:N_c
+            push!(eqs, n_node[j] * xᵢ_node[j, i] - n_node[j - 1] * xᵢ_node[j - 1, i] ~
+                       R_node[j, i] * V_seg)
+        end
+
+        # Total molar flow (sum of component flows)
+        push!(eqs, n_node[j] ~ sum([n_node[j] * xᵢ_node[j, i] for i in 1:N_c]))
+
+        # Pressure drop (simplified: assume small or specify separately)
+        # For now, assume isobaric
+        push!(eqs, p_node[j] ~ p_node[j - 1])
+
+        # Energy balance: d(n*h)/dz = ΔHᵣ * A + Q
+        push!(eqs, n_node[j] * h_node[j] - n_node[j - 1] * h_node[j - 1] ~
+                   ΔHᵣ_node[j] + Q_node[j])
+    end
+
+    # Normalize mole fractions (ensure sum = 1)
+    for j in 2:N_nodes
+        push!(eqs, 1.0 ~ sum([xᵢ_node[j, i] for i in 1:N_c]))
+    end
+
+    # Outlet conditions
+    push!(eqs, F_out ~ n_node[N_nodes])
+    push!(eqs, T_out ~ T_node[N_nodes])
+
+    # Outlet connector
+    push!(eqs, Out.n ~ -F_out)
+    push!(eqs, Out.T ~ T_out)
+    push!(eqs, Out.p ~ p_node[N_nodes])
+    push!(eqs, Out.ϱ ~ ϱ_node[N_nodes])
+    push!(eqs, Out.h ~ h_node[N_nodes])
+    for i in 1:N_c
+        push!(eqs, Out.xᵢ[i] ~ xᵢ_node[N_nodes, i])
+    end
+
+    return ODESystem(eqs, t, collect(Iterators.flatten(vars)),
+        collect(Iterators.flatten(pars));
+        name, systems = [In, Out])
+end
+
+"""
+    SteadyStatePFR(ms::MaterialSource; name, L, A, N_nodes, i_reacts, adiabatic)
+
+Create a steady-state Plug Flow Reactor. This is a convenience wrapper that
+sets up the PFR for steady-state analysis (d/dt = 0).
+
+All arguments are the same as PFR.
+"""
+@component function SteadyStatePFR(ms::MaterialSource;
+        name,
+        L::Float64,
+        A::Float64,
+        N_nodes::Int = 20,
+        i_reacts = 1:length(ms.reaction),
+        adiabatic::Bool = true)
+    # For steady state, we simply create the PFR - the time derivatives
+    # will be zero in the discretized spatial equations
+    return PFR(ms; name, L, A, N_nodes, i_reacts, adiabatic)
+end
+
+"""
+    SimplePFR(ms::MaterialSource; name, L, A, Q_total)
+
+A simplified PFR model without spatial discretization. Uses a single
+control volume with residence time based on volume and flow rate.
+Suitable for preliminary calculations.
+
+# Arguments
+- `ms::MaterialSource`: Material source
+- `name`: System name
+- `L::Float64`: Reactor length (m)
+- `A::Float64`: Cross-sectional area (m²)
+- `Q_total::Float64`: Total volumetric flow rate (m³/s) - optional, can be from inlet
+
+# Returns
+An `ODESystem` for the simplified PFR.
+"""
+@component function SimplePFR(ms::MaterialSource;
+        name,
+        L::Float64,
+        A::Float64,
+        i_reacts = 1:length(ms.reaction))
+    N_c = ms.N_c
+    N_r = length(ms.reaction[i_reacts])
+    V_reactor = L * A
+
+    # Connectors
+    @named In = MaterialConnector(ms)
+    @named Out = MaterialConnector(ms)
+
+    # Parameters
+    pars = @parameters begin
+        V = V_reactor, [description = "Reactor volume (m³)"]
+        L_reactor = L, [description = "Reactor length (m)"]
+    end
+
+    # Variables
+    vars = @variables begin
+        τ(t), [description = "Residence time (s)"]
+        T(t), [description = "Outlet temperature (K)"]
+        p(t), [description = "Pressure (Pa)"]
+        (xᵢ(t))[1:N_c], [description = "Outlet mole fractions"]
+        ϱ(t), [description = "Molar density (mol/m³)"]
+        h(t), [description = "Molar enthalpy (J/mol)"]
+        n_out(t), [description = "Outlet molar flow (mol/s)"]
+        (r(t))[1:N_r], [description = "Reaction rates (mol/s/m³)"]
+        (R(t))[1:N_c], [description = "Species production rates (mol/s/m³)"]
+        ΔHᵣ(t), [description = "Heat of reaction (J/s)"]
+    end
+
+    eqs = Equation[]
+
+    # Residence time
+    push!(eqs, τ ~ V * ϱ / In.n)
+
+    # Pressure (isobaric)
+    push!(eqs, p ~ In.p)
+
+    # Density
+    push!(eqs, ϱ ~ ms.molar_density(p, T, collect(xᵢ); phase = "unknown"))
+
+    # Enthalpy
+    push!(eqs, h ~ ms.VT_enthalpy(ϱ, T, collect(xᵢ)))
+
+    # Reaction rates (using average composition)
+    for (ir, reac) in enumerate(ms.reaction[i_reacts])
+        # Use inlet composition for rate calculation (first-order approximation)
+        push!(eqs, r[ir] ~ reac.r(p, T, collect(In.xᵢ)))
+    end
+
+    # Species production rates
+    for i in 1:N_c
+        push!(eqs, R[i] ~ sum([r[ir] * ms.reaction[i_reacts][ir].ν[i] for ir in 1:N_r]))
+    end
+
+    # Component balance (steady state): F_in * x_in + R * V = F_out * x_out
+    for i in 1:N_c
+        push!(eqs, In.n * In.xᵢ[i] + R[i] * V ~ n_out * xᵢ[i])
+    end
+
+    # Total molar balance
+    push!(eqs, n_out ~ In.n + sum([R[i] * V for i in 1:N_c]))
+
+    # Heat of reaction
+    push!(eqs, ΔHᵣ ~ sum([r[ir] * ms.reaction[i_reacts][ir].Δhᵣ(T) * V for ir in 1:N_r]))
+
+    # Energy balance (adiabatic): F_in * h_in = F_out * h_out + ΔHᵣ
+    push!(eqs, In.n * In.h ~ n_out * h + ΔHᵣ)
+
+    # Mole fraction normalization
+    push!(eqs, 1.0 ~ sum([xᵢ[i] for i in 1:N_c]))
+
+    # Outlet connector
+    push!(eqs, Out.n ~ -n_out)
+    push!(eqs, Out.T ~ T)
+    push!(eqs, Out.p ~ p)
+    push!(eqs, Out.ϱ ~ ϱ)
+    push!(eqs, Out.h ~ h)
+    for i in 1:N_c
+        push!(eqs, Out.xᵢ[i] ~ xᵢ[i])
+    end
+
+    return ODESystem(eqs, t, collect(Iterators.flatten(vars)),
+        collect(Iterators.flatten(pars));
+        name, systems = [In, Out])
+end
+
+export PFR, SteadyStatePFR, SimplePFR
diff --git a/src/reactors/reaction_sets.jl b/src/reactors/reaction_sets.jl
new file mode 100644
index 0000000..8f66bd7
--- /dev/null
+++ b/src/reactors/reaction_sets.jl
@@ -0,0 +1,195 @@
+"""
+    Reaction sets and kinetics utilities for ProcessSimulator.jl
+
+This module provides data structures and utilities for managing sets of
+chemical reactions and computing reaction kinetics.
+"""
+
+"""
+    ReactionSet
+
+A collection of related reactions that occur together in a reactor.
+Provides convenient access to stoichiometry, kinetics, and thermodynamic
+properties for multiple reactions.
+
+# Fields
+- `reactions::Vector{Reaction}`: Vector of individual reactions
+- `name::String`: Name/identifier for the reaction set
+- `components::Vector{String}`: Names of all components involved
+"""
+struct ReactionSet
+    reactions::Vector{Reaction}
+    name::String
+    components::Vector{String}
+
+    function ReactionSet(reactions::Vector{Reaction}, components::Vector{String};
+                         name::String = "ReactionSet")
+        # Validate that all reactions have consistent stoichiometry length
+        N_c = length(components)
+        for (i, r) in enumerate(reactions)
+            if length(r.ν) != N_c
+                throw(ArgumentError(
+                    "Reaction $i has stoichiometry of length $(length(r.ν)), " *
+                    "expected $N_c (number of components)"))
+            end
+        end
+        new(reactions, name, components)
+    end
+end
+
+"""
+    ReactionSet(; reactions, components, name="ReactionSet")
+
+Keyword constructor for ReactionSet.
+"""
+function ReactionSet(; reactions::Vector{Reaction}, components::Vector{String},
+                     name::String = "ReactionSet")
+    ReactionSet(reactions, components; name = name)
+end
+
+"""
+    nreactions(rs::ReactionSet)
+
+Return the number of reactions in the set.
+"""
+nreactions(rs::ReactionSet) = length(rs.reactions)
+
+"""
+    ncomponents(rs::ReactionSet)
+
+Return the number of components in the reaction set.
+"""
+ncomponents(rs::ReactionSet) = length(rs.components)
+
+"""
+    stoichiometry_matrix(rs::ReactionSet)
+
+Return the stoichiometric matrix (N_reactions × N_components).
+Each row corresponds to a reaction, each column to a component.
+Negative values indicate reactants, positive values indicate products.
+"""
+function stoichiometry_matrix(rs::ReactionSet)
+    N_r = nreactions(rs)
+    N_c = ncomponents(rs)
+    ν = zeros(N_r, N_c)
+    for (i, r) in enumerate(rs.reactions)
+        ν[i, :] .= r.ν
+    end
+    return ν
+end
+
+"""
+    reaction_rates(rs::ReactionSet, p, T, x)
+
+Compute the reaction rates for all reactions in the set.
+Returns a vector of rates in mol/s.
+
+# Arguments
+- `rs::ReactionSet`: The reaction set
+- `p`: Pressure (Pa)
+- `T`: Temperature (K)
+- `x`: Mole fractions (vector of length N_components)
+"""
+function reaction_rates(rs::ReactionSet, p, T, x)
+    return [r.r(p, T, x) for r in rs.reactions]
+end
+
+"""
+    species_production_rates(rs::ReactionSet, p, T, x)
+
+Compute the net production rate of each species due to all reactions.
+Returns a vector of rates in mol/s for each component.
+
+# Arguments
+- `rs::ReactionSet`: The reaction set
+- `p`: Pressure (Pa)
+- `T`: Temperature (K)
+- `x`: Mole fractions (vector of length N_components)
+"""
+function species_production_rates(rs::ReactionSet, p, T, x)
+    rates = reaction_rates(rs, p, T, x)
+    ν = stoichiometry_matrix(rs)
+    return ν' * rates  # N_c × 1 vector
+end
+
+"""
+    reaction_enthalpies(rs::ReactionSet, T)
+
+Compute the enthalpy of reaction for each reaction at temperature T.
+Returns a vector of enthalpies in J/mol.
+"""
+function reaction_enthalpies(rs::ReactionSet, T)
+    return [r.Δhᵣ(T) for r in rs.reactions]
+end
+
+"""
+    total_heat_of_reaction(rs::ReactionSet, p, T, x)
+
+Compute the total heat released/absorbed by all reactions.
+Returns heat rate in J/s.
+
+# Arguments
+- `rs::ReactionSet`: The reaction set
+- `p`: Pressure (Pa)
+- `T`: Temperature (K)
+- `x`: Mole fractions (vector of length N_components)
+"""
+function total_heat_of_reaction(rs::ReactionSet, p, T, x)
+    rates = reaction_rates(rs, p, T, x)
+    enthalpies = reaction_enthalpies(rs, T)
+    return sum(rates .* enthalpies)
+end
+
+"""
+    ArrheniusKinetics
+
+Helper struct for Arrhenius-type reaction kinetics:
+    k = A * exp(-Ea / (R * T))
+
+# Fields
+- `A`: Pre-exponential factor (units depend on reaction order)
+- `Ea`: Activation energy (J/mol)
+"""
+struct ArrheniusKinetics
+    A::Float64   # Pre-exponential factor
+    Ea::Float64  # Activation energy (J/mol)
+end
+
+"""
+    rate_constant(k::ArrheniusKinetics, T)
+
+Compute the rate constant at temperature T (K) using Arrhenius equation.
+"""
+function rate_constant(k::ArrheniusKinetics, T)
+    R = 8.314  # J/(mol·K)
+    return k.A * exp(-k.Ea / (R * T))
+end
+
+"""
+    PowerLawKinetics
+
+Helper struct for power-law rate expressions:
+    r = k(T) * prod(C_i^order_i)
+
+# Fields
+- `arrhenius::ArrheniusKinetics`: Temperature-dependent rate constant
+- `orders::Vector{Float64}`: Reaction orders for each component
+"""
+struct PowerLawKinetics
+    arrhenius::ArrheniusKinetics
+    orders::Vector{Float64}
+end
+
+"""
+    reaction_rate(k::PowerLawKinetics, T, C)
+
+Compute the reaction rate given temperature T and concentrations C.
+"""
+function reaction_rate(k::PowerLawKinetics, T, C)
+    k_val = rate_constant(k.arrhenius, T)
+    return k_val * prod(C .^ k.orders)
+end
+
+export ReactionSet, nreactions, ncomponents, stoichiometry_matrix
+export reaction_rates, species_production_rates, reaction_enthalpies, total_heat_of_reaction
+export ArrheniusKinetics, PowerLawKinetics, rate_constant, reaction_rate
diff --git a/src/reactors/types.jl b/src/reactors/types.jl
new file mode 100644
index 0000000..9719d11
--- /dev/null
+++ b/src/reactors/types.jl
@@ -0,0 +1,46 @@
+"""
+    Abstract type hierarchy for reactors in ProcessSimulator.jl
+
+The reactor type hierarchy provides a structured way to categorize different
+reactor models based on their fundamental modeling approach:
+
+- `AbstractReactor`: Base type for all reactor models
+- `AbstractEquilibriumReactor`: Reactors based on thermodynamic equilibrium (Gibbs minimization)
+- `AbstractKineticReactor`: Reactors based on reaction kinetics (rate equations)
+"""
+
+"""
+    AbstractReactor
+
+Base abstract type for all reactor components. Reactors transform input streams
+through chemical reactions to produce output streams.
+"""
+abstract type AbstractReactor end
+
+"""
+    AbstractEquilibriumReactor <: AbstractReactor
+
+Abstract type for equilibrium-based reactors. These reactors assume chemical
+equilibrium is reached and determine the outlet composition by minimizing
+Gibbs free energy subject to elemental balance constraints.
+
+Examples include:
+- GibbsReactor: General equilibrium reactor using Gibbs minimization
+- EquilibriumReactor: Simplified equilibrium reactor with specified reactions
+"""
+abstract type AbstractEquilibriumReactor <: AbstractReactor end
+
+"""
+    AbstractKineticReactor <: AbstractReactor
+
+Abstract type for kinetics-based reactors. These reactors model the rate
+of chemical reactions using rate expressions (e.g., Arrhenius kinetics).
+
+Examples include:
+- CSTR (Continuous Stirred Tank Reactor): Perfect mixing assumption
+- PFR (Plug Flow Reactor): No axial mixing, 1D spatial variation
+- BatchReactor: Time-varying composition in closed vessel
+"""
+abstract type AbstractKineticReactor <: AbstractReactor end
+
+export AbstractReactor, AbstractEquilibriumReactor, AbstractKineticReactor
diff --git a/test/reactors/reactor_types_test.jl b/test/reactors/reactor_types_test.jl
new file mode 100644
index 0000000..e7b6904
--- /dev/null
+++ b/test/reactors/reactor_types_test.jl
@@ -0,0 +1,50 @@
+using ProcessSimulator
+using Test
+
+const PS = ProcessSimulator
+
+@testset "Reactor Type Hierarchy" begin
+    # Test abstract type hierarchy
+    @test PS.AbstractEquilibriumReactor <: PS.AbstractReactor
+    @test PS.AbstractKineticReactor <: PS.AbstractReactor
+end
+
+@testset "ReactionSet Utilities" begin
+    # Define a simple reaction: A -> B
+    r1 = PS.Reaction(
+        ν = [-1.0, 1.0],
+        r = (p, T, x) -> 1e-3 * x[1],  # First order in A
+        Δhᵣ = (T) -> -50000.0  # Exothermic, 50 kJ/mol
+    )
+
+    rs = PS.ReactionSet([r1], ["A", "B"]; name = "SimpleReaction")
+
+    @test PS.nreactions(rs) == 1
+    @test PS.ncomponents(rs) == 2
+
+    ν_matrix = PS.stoichiometry_matrix(rs)
+    @test ν_matrix[1, 1] == -1.0  # A is consumed
+    @test ν_matrix[1, 2] == 1.0   # B is produced
+end
+
+@testset "ArrheniusKinetics" begin
+    # Create Arrhenius kinetics for a typical reaction
+    # A = 1e10 s^-1, Ea = 50 kJ/mol
+    k = PS.ArrheniusKinetics(1e10, 50000.0)
+
+    # Test rate constant at different temperatures
+    k_300 = PS.rate_constant(k, 300.0)
+    k_400 = PS.rate_constant(k, 400.0)
+
+    # Higher temperature should give higher rate constant
+    @test k_400 > k_300
+
+    # k = A * exp(-Ea/(R*T)) where R = 8.314 J/(mol·K)
+    # At 300K: k = 1e10 * exp(-50000/(8.314*300)) ≈ 1e10 * exp(-20.05) ≈ 19.7
+    expected_k_300 = 1e10 * exp(-50000.0 / (8.314 * 300.0))
+    @test isapprox(k_300, expected_k_300, rtol = 1e-6)
+
+    # At 400K: k = 1e10 * exp(-50000/(8.314*400)) ≈ 1e10 * exp(-15.03) ≈ 2.97e3
+    expected_k_400 = 1e10 * exp(-50000.0 / (8.314 * 400.0))
+    @test isapprox(k_400, expected_k_400, rtol = 1e-6)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 437fc5a..3e555aa 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,3 +9,7 @@ end
 @safetestset "Reactors" begin
     include("reactors/simple_cstr.jl")
 end
+
+@safetestset "Reactor Types" begin
+    include("reactors/reactor_types_test.jl")
+end