Hello,
I am trying to link the internal heat source to the reaction rate of pyrite to simulate an exothermic reaction. I created a 64 bit library and was able load the function when into X1t. When the simulation starts, the program instantly closes. Any suggestions about how I broke this?
Thanks,
Peter
#X1t input file
data = thermo.tdat verify
conductivity = conductivity-USGS.dat
suppress ALL
unsuppress Hematite K-feldspar Pyrite
interval start at 0 day, fluid = fluid_1
time end at .25 m.y.
length = 250 m
width = 1 m
height = 1 m
Nx = 25
discharge start = .307 m3/m2/yr
diffusion_coef = 1e-7 cm2/s
temperature = 100 C
heat_source = file D:\Model\Teena\Heat\heat.dll:heat J/cm3/yr transient, Tmax = 300 C
scope = initial
H2O = 1 free kg
Na+ = 1 g/kg
balance on Cl-
Fe++ = 1 mg/kg
swap Pyrite for O2(aq)
Pyrite = 10 free volume%
SO4-- = 20 mg/kg
pH = 6.5
scope = fluid_1
H2O = 1 free kg
Na+ = 1 g/kg
balance on Cl-
Fe++ = 1 mg/kg
swap e- for O2(aq)
Eh = .2 V
SO4-- = 20 mg/kg
pH = 6
scope
kinetic Pyrite pre-exp = 3.02e-12 act_eng = 56900 surface = 1e4
// Heatlib.cpp : Defines the functions for the library.
//
#define EXPORT extern "C" __declspec(dllexport)
#define NOMINMAX
#include <stdio.h>
#include <math.h>
#include "gwb_context.h"
// heat: Examine the amount of pyrite oxidized and return a heat source
// proportional to the pyrite reaction rate and the time step.
EXPORT double heat(int i, Nodal_block &n, GWBcontext &c)
{
double pre, act_eng, qnk, surf, temp;
pre = n.rct[0].preexp;
act_eng = n.rct[0].acteng;
qnk = n.rct[0].QoverK;
surf = n.rct[0].surfac;
temp = n.tempc + 273.15;
if (qnk > 1.0)
return 0.0;
//This is the reaction rate in mol/s which needs to be in J/cm3/s
//so divide by the domain size and number of nodal blocks
//which gets mol/cm3/s and then multiply by the energy of the reaction
//in J/mol which I took from https://doi.org/10.3390/min11060565
return pre*surf*(1.0-qnk)*exp(-act_eng/8.314/temp)/c.Length/c.Width/c.Height/c.Nnode*1409000.0;
}