Material coefficients with time derivatives of field variables
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', }
Hi Max,
On 16/06/2021 23:37, Max Maximilian wrote:
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.
By that you mean that liquid_saturation_dt stays zero?
There might be a bug that manifests in multiphysics problems (with several variables and time derivatives). I will have a look.
Does it mean I need to reformulate the equations that exclude time derivative terms?
If you could avoid using variables with history, that might help in the meantime.
Thanks for the report! r.
Hi Robert,
that's it. liquid_saturation_dt and rho_dt stay zero.
I will try reducing variables with history, then update a feedback.
Thanks a lot!
Hey Robert,
I've got trouble for field-variable-dependent matrix.
For example, if a stiffness matrix is strongly dependent on u and grad u, then I use the way shown in examples, call pb.evaluate(). It seems that it cannot work well when the pdes are too complex (multiphase multiple component flows).
Are there any other ways to calculate variable-/grad_variable dependent matrix?
Hi Max,
On 22/06/2021 22:22, Max Maximilian wrote:
Hey Robert,
I've got trouble for field-variable-dependent matrix.
For example, if a stiffness matrix is strongly dependent on u and grad u, then I use the way shown in examples, call pb.evaluate(). It seems that it cannot work well when the pdes are too complex (multiphase multiple component flows).
By "cannot work well" you mean that the iterations do not converge?
Are there any other ways to calculate variable-/grad_variable dependent matrix?
There are terms for nonlinear elasticity (hyperelasticity - see e.g. dw_tl_* terms). In a similar manner, other nonlinear relations can be differentiated (analytically by hand/sympy/... - there is no automatic differentiation code included) to have consistent tangents, and then use the Newton's solver. This, however, can also have issues with convergence. Usually some step reduction/line search is required.
r.
participants (3)
-
Mary Regner -
Max Maximilian -
Robert Cimrman