Hey all,
I found my script cannot successfully compute the volume_dot matrix including time derivatives of field variables, since it stayed 0 at all time.
Does it mean I need to reformulate the equations that exclude time derivative terms?
## CODE
from __future__ import absolute_import
from sfepy.examples.dg.example_dg_common import *
from sfepy import data_dir
from sfepy.base.base import output
import numpy as nm
from sfepy.discrete.fem import Mesh
from sfepy.discrete.fem.meshio import UserMeshIO
import matplotlib.pyplot as plt
from sfepy.base.base import Struct
filename_mesh=None
X1 = 0
XN = 2e-4
n_nod = 2001
n_el = n_nod - 1
if filename_mesh is None:
filename_mesh = get_gen_1D_mesh_hook(X1, XN, n_nod)
approx_order=2
t0=0.0
t1=0.5
porosity = 0.52 # PTL porosity
theta_PTL = nm.pi*30./180. # PTL contact angle
KG = 4.5e-15 # Permeability, m2
R = 8.3145 # Ideal gas constant, J/mol/K
T = 273.15+80 # Ambient temperature, K
patm = 101325.0 # Ambient pressure
sCh = 0.5 # Channel liquid saturation
sl = 0.4 # Left boundary liquid saturation
Mox = 0.032 # Molar mass density of oxygen, kg/mol
i_app = 12000. # Applied current density, A/m2
F = 96485 # Faraday constant, C/mol
Mw = 0.018 # Molar mass density of water, kg/mol
muG = 1.881e-5 # Gas viscosity, Pa*s
muL = 3.56e-4 # Liqiud viscosity, Pa*s
rho_H2O = 996 # Water density, kg/m3
m = 3 # Relative permeabiltiy coefficient
psat = 10**(8.07131-1730.63/(233.426+T-273.15))*133.322 # Saturated water gas pressure, Pa
kc = 1e4 # Condensation rate coefficient, 1/s
kv = 5.1e-5 # Evaporation rate coefficient, 1/(Pa*s)
Dox = 2.1e-5 # Diffusivity of oxygen, m^2/s
gas_rate = i_app/4/F*Mox
liquid_rate = -i_app/2/F*Mw
# Boundary & Initial pressure/mass fraction
pc0 = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(1.417*(1-sCh)-2.120*(1-sCh)**2+1.263*(1-sCh)**3)
rho00 = Mox*(patm+pc0-psat)/R/T+psat*Mw/R/T
c0 = Mox*(patm+pc0-psat)/(Mox*(patm+pc0-psat)+psat*Mw)
pcl = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(1.417*(1-0.4)-2.120*(1-0.4)**2+1.263*(1-0.4)**3)
rhol = Mox*(patm+pcl-psat)/R/T+psat*Mw/R/T
cl = Mox*(patm+pcl-psat)/(Mox*(patm+pcl-psat)+psat*Mw)
example_name = "fluid_dynamics_1D"
def get_ic_w(coors, ic=None):
pc = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(1.417*(1-((sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl))-2.120*(1-((sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl))**2+1.263*(1-((sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl))**3)
rho0 = Mox*(patm+pc-psat)/R/T+psat*Mw/R/T
rho0.shape = (coors.shape[0], 1, 1)
return rho0
def get_ic_s(coors, ic=None):
s0 = (sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl
s0.shape = (coors.shape[0], 1, 1)
return s0
def get_ic_c(coors, ic=None):
pc = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(1.417*(1-((sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl))-2.120*(1-((sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl))**2+1.263*(1-((sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl))**3)
c0 = (Mox*(patm+pc-psat))/(Mox*(patm+pc-psat)+psat*Mw)
c0.shape = (coors.shape[0], 1, 1)
return c0
def post_process(out, pb, state, extend=False):
# Field variable evalution
rho_grad = pb.evaluate('ev_grad.i.Omega(rho)',
mode='el_avg', verbose=False)
s_grad = pb.evaluate('ev_grad.i.Omega(s)',
mode='el_avg', verbose=False)
c_grad = pb.evaluate('ev_grad.i.Omega(c)',
mode='el_avg', verbose=False)
liquid_saturation_values = pb.evaluate('ev_volume_integrate.i.Omega(s)',
mode='el_avg', verbose=False)
rho_values = pb.evaluate('ev_volume_integrate.i.Omega(rho)',
mode='el_avg', verbose=False)
conc_values = pb.evaluate('ev_volume_integrate.i.Omega(c)',
mode='el_avg', verbose=False)
dpcds = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(-1.417+2.120*2*(1-liquid_saturation_values)-3*1.263*(1-liquid_saturation_values)**2)
liquid_velocity = -KG*liquid_saturation_values**m/muL*R*T*((conc_values/Mox+(1-conc_values)/Mw)*rho_grad+(1/Mox-1/Mw)*rho_values*c_grad-dpcds*s_grad/R/T)
gas_velocity = -KG*(1-liquid_saturation_values)**m/muG*R*T*((conc_values/Mox+(1-conc_values)/Mw)*rho_grad+(1/Mox-1/Mw)*rho_values*c_grad)
out['liquid_velocity'] = Struct(name='output_data',
mode='cell', data=liquid_velocity,
dofs=None)
out['gas_velocity'] = Struct(name='output_data',
mode='cell', data=gas_velocity,
dofs=None)
return out
def get_coef(ts, coors, problem, equations=None, mode=None, **kwargs):
"""
Calculates the coefficient matrix and returns them.
This relation results in velocities and diffusivity.
"""
if mode == 'qp':
# Field and gradient values in quadrature points coordinates given by integral i
# - they are the same as in `coors` argument.
if ts.step == 0:
liquid_saturation_values = (sCh-sl)/(XN-X1)*(coors[:,0]-X1)+sl
else:
liquid_saturation_values = problem.evaluate('ev_volume_integrate.i.Omega(s)',
mode='qp', verbose=False)
liquid_saturation_values.shape = (coors.shape[0], 1, 1)
if ts.step == 0:
liquid_saturation_dt = nm.zeros(coors.shape[0])
else:
liquid_saturation_dt = problem.evaluate('ev_volume_integrate.i.Omega(ds/dt)',
mode='qp', verbose=False)
liquid_saturation_dt.shape = (coors.shape[0], 1, 1)
if ts.step == 0:
dpcds = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(-1.417+2.120*2*(1-liquid_saturation_values)-3*1.263*(1-liquid_saturation_values)**2)
grad_rho_values = Mox/R/T*dpcds*(sCh-sl)/(XN-X1)
else:
grad_rho_values = problem.evaluate('ev_grad.i.Omega(rho)',
mode='qp', verbose=False)
grad_rho_values.shape = (coors.shape[0], 1, 1)
if ts.step == 0:
pc = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(1.417*(1-liquid_saturation_values)-2.120*(1-liquid_saturation_values)**2+1.263*(1-liquid_saturation_values)**3)
rho_values = Mox*(patm+pc-psat)/R/T+psat*Mw/R/T
else:
rho_values = problem.evaluate('ev_volume_integrate.i.Omega(rho)',
mode='qp', verbose=False)
rho_values.shape = (coors.shape[0], 1, 1)
if ts.step == 0:
rho_dt = nm.zeros(coors.shape[0])
else:
rho_dt = problem.evaluate('ev_volume_integrate.i.Omega(drho/dt)',
mode='qp', verbose=False)
rho_dt.shape = (coors.shape[0], 1, 1)
if ts.step == 0:
conc_values = (Mox*(patm+pc-psat))/(Mox*(patm+pc-psat)+psat*Mw)
else:
conc_values = problem.evaluate('ev_volume_integrate.i.Omega(c)',
mode='qp', verbose=False)
conc_values.shape = (coors.shape[0], 1, 1)
if ts.step == 0:
grad_conc_values = psat*Mw/(Mox*(patm+pc-psat)+psat*Mw)**2*Mox*dpcds*(sCh-sl)/(XN-X1)
else:
grad_conc_values = problem.evaluate('ev_grad.i.Omega(c)',
mode='qp', verbose=False)
grad_conc_values.shape = (coors.shape[0], 1, 1)
# Gas velocity
u1 = -KG*(1-liquid_saturation_values)**m/muG*R*T*((conc_values/Mox+(1-conc_values)/Mw)*grad_rho_values+(1/Mox-1/Mw)*rho_values*grad_conc_values)
## Gas density mass conservation
# Mass matrix of gas density: (1-s)*eps
mass_dens = (1-liquid_saturation_values)*porosity
# Source matrix of gas density: porosity*ds/dt
source_dens = liquid_saturation_dt*porosity
# Diffusion coupling matrix of gas density: u1
diffcop_dens = u1
# Reaction matrix of gas density: Rsat (Phase transistion)
p = rho_values*(1-conc_values)*R*T/Mw+conc_values*rho_values*R*T/Mox
reac_dens = kc*porosity*(1-liquid_saturation_values)*Mw/R/T*(rho_values*(1-conc_values)*R*T/Mw-psat)*nm.heaviside(rho_values*(1-conc_values)*R*T/Mw-psat, 0) + kv*porosity*liquid_saturation_values*rho_H2O*(rho_values*(1-conc_values)*R*T/Mw-psat)*nm.heaviside(psat-rho_values*(1-conc_values)*R*T/Mw, 0)
## Liquid water mass conservation
# Mass matrix of liquid saturation: porosity*rho_H2O
mass_sat = rho_H2O*porosity*nm.ones(coors.shape[0])
mass_sat.shape = (coors.shape[0], 1, 1)
# Diffusion_r matrix of liquid saturation: rho_H2O*K*s^m/muL*R*T*((c/Mox+(1-c)/Mw)*grad_rho+grad_c*(1/Mox-1/Mw)*rho)
Diff_r_sat = rho_H2O*KG*liquid_saturation_values**m/muL*R*T*((conc_values/Mox+(1-conc_values)/Mw)*grad_rho_values+(1/Mox-1/Mw)*grad_conc_values*rho_values)
# Stiffness matrix of liquid saturation: rho_H2O*dpcds*K*s^m/muL
dpcds = 235.8e-3*(1-T/647.096)**1.256*(1-0.625*(1-T/647.096))*nm.cos(theta_PTL)*nm.sqrt(porosity/KG)*(-1.417+2.120*2*(1-liquid_saturation_values)-3*1.263*(1-liquid_saturation_values)**2)
stiff_sat = rho_H2O*KG*liquid_saturation_values**m/muL*dpcds
# Reaction matrix of liquid saturation: Rsat
reac_sat = reac_dens
## Mass contration of Oxygen
# Mass matrix of oxygen concentration
mass_conc = porosity*(1-liquid_saturation_values)*rho_values
# Source 1 matrix of oxygen concentration: eps*(1-s)*drho/dt
source1_conc = rho_dt*porosity*(1-liquid_saturation_values)
# Source 2 matrix of oxygen concentration: eps*ds/dt*rho
source2_conc = liquid_saturation_dt*porosity*rho_values
# Diffusion coupling matrix of oxygen concentration: u1*rho
diffcop_conc = u1*rho_values
# Stiffness matrix of oxygen concentration: rho*Deff
stiff_conc = rho_values*Dox*(porosity*(1-liquid_saturation_values))**1.5
output('min:', source_dens.min(), 'max:', source_dens.max())
return {'mass_dens' : mass_dens, 'source_dens' : source_dens, 'diffcop_dens' : diffcop_dens, 'reac_dens' : reac_dens,
'mass_sat' : mass_sat, 'Diff_r_sat' : Diff_r_sat, 'stiff_sat' : stiff_sat, 'reac_sat' : reac_sat,
'mass_conc' : mass_conc, 'source1_conc' : source1_conc, 'source2_conc' : source2_conc, 'diffcop_conc' : diffcop_conc, 'stiff_conc' : stiff_conc}
materials = {
'coef' : 'get_coef',
'flux' : ({
'gas_rate' : gas_rate,
'liquid_rate' : liquid_rate,
'conc_rate' : 0,
},),
}
fields = {
'gas_density': ('real', 'scalar', 'Omega', approx_order),
'liquid_saturation': ('real', 'scalar', 'Omega', approx_order),
'gas_mass_fraction': ('real', 'scalar', 'Omega', approx_order),
}
variables = {
'rho': ('unknown field', 'gas_density', 0, 1),
'v1': ('test field', 'gas_density', 'rho'),
's': ('unknown field', 'liquid_saturation', 1, 1),
'v2': ('test field', 'liquid_saturation', 's'),
'c': ('unknown field', 'gas_mass_fraction', 2, 1),
'v3': ('test field', 'gas_mass_fraction', 'c'),
}
regions = {
'Omega' : 'all',
'Gamma_Left' : ('vertices in (x < 1e-7)', 'facet'),
'Gamma_Right' : ('vertices in (x > 1.999e-4)', 'facet'),
}
ebcs = {
'right' : ('Gamma_Right', {'rho.0' : rho00, 's.0' : sCh, 'c.0' : c0}),
'left' : ('Gamma_Left', {'rho.0' : rhol, 's.0' : 0.4, 'c.0' : cl}),
}
integrals = {
'i' : 2*approx_order,
}
equations = {
'gas_density' : """dw_volume_dot.i.Omega(coef.mass_dens, v1, drho/dt) - dw_volume_dot.i.Omega(coef.source_dens, v1, rho) -
dw_diffusion_coupling.i.Omega(coef.diffcop_dens, v1, rho) + dw_volume_integrate.i.Omega(coef.reac_dens, v1) = 0""",
'liquid_saturation' : """dw_volume_dot.i.Omega(coef.mass_sat, v2, ds/dt) + dw_diffusion_r.i.Omega(coef.Diff_r_sat, v2)
- dw_laplace.i.Omega(coef.stiff_sat, v2, s) - dw_volume_integrate.i.Omega(coef.reac_sat, v2) = 0""",
'gas_mass_fraction' : """dw_volume_dot.i.Omega(coef.mass_conc, v3, dc/dt) + dw_volume_dot.i.Omega(coef.source1_conc, v3, c) - dw_volume_dot.i.Omega(coef.source2_conc, v3, c)
- dw_diffusion_coupling.i.Omega(coef.diffcop_conc, v3, c) + dw_laplace.i.Omega(coef.stiff_conc, v3, c) = 0 """,
}
# -dw_surface_ndot.i.Gamma_Left(flux.gas_rate, v1) -dw_surface_ndot.i.Gamma_Left(flux.liquid_rate, v2) -dw_surface_ndot.i.Gamma_Left(flux.conc_rate, v3)
functions = {
'get_coef' : (get_coef,),
'get_ic_w' : (get_ic_w,),
'get_ic_s' : (get_ic_s,),
'get_ic_c' : (get_ic_c,),
}
ics = {
'ic' : ('Omega', {'rho.0' : 'get_ic_w', 's.0' : 'get_ic_s', 'c.0' : 'get_ic_c'}),
}
solvers = {
'ls' : ('ls.mumps', {
}),
'newton' : ('nls.newton', {
'eps_a' : 1e-4,
'eps_r' : 1e-3,
'i_max' : 5,
}),
'tss' : ('ts.adaptive', {
't0' : t0,
't1' : t1,
'dt' : None,
'verbose' : 1,
'quasistatic' : False,
'n_step' : 100,
}),
}
options = {
'ts' : 'tss',
'nls' : 'newton',
'ls' : 'ls',
'save_times' : 10,
'output_dir' : 'output/' + example_name,
'output_format' : "vtk",
'post_process_hook' : 'post_process',
}
