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 outdef 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', }