SfePy

sfepy@python.org

6 participants
839 discussions

Importing Stress-strain behavior rather than solving PDE

by kshargh.ali＠gmail.com

Hello all, It's only a while since I started to work with SfePy and It's awesome. I also have a concern and I really appreciate it if you could please share your thoughts on it with me: What I want to do is to import stress-strain curve as a table ( or function ) and do the tensile test on a sample rather than defining Young's modulus and Poisson's ratio to the software and let it solve constitutive equation to obtain the result of tensile test. Is it something I can do with SfePy ? Best, Ali

1 day, 10 hours

Interactive Multi-Material Meshes

by Joshua colletta

Hi all, I'm trying to model a magnet encased in a cube of elastomer for part of a project. I've seemingly tried everything under the sun, but ive been rather limited due to the poor interactive documentation. I have to core problems: Defining a region by a function: Material parameters can be defined over just a region, as such a function to define an inner regions of given bounds: def inner_region_(coors, domain=None ,magnet=None): x, y, z = coors[:, 0], coors[:, 1], coors[:, 2] mag_min_x, mag_max_x = str((magnet[1][0] - (magnet[0][0] / 2))), str((magnet[1][0] + (magnet[0][0] / 2))) mag_min_y, mag_max_y = str((magnet[1][1] - (magnet[0][1] / 2))), str((magnet[1][1] + (magnet[0][1] / 2))) mag_min_z, mag_max_z = str((magnet[1][2] - (magnet[0][2] / 2))), str((magnet[1][2] + (magnet[0][2] / 2))) parse_def_min = 'vertices in (' + x + ' > ' + mag_min_x + ') & (' + y + ' > ' + mag_min_y + ') & (' + z + ' > ' + mag_min_z + ')' parse_def_max = '& (' + x + ' < ' + mag_max_x + ') & (' + y + ' < ' + mag_max_y + ') & (' + z + ' < ' + mag_max_z + ')' #pls dont hate my ugly strings :) parse_def = parse_def_min + parse_def_max flag = nm.where(parse_def)[0] print(flag) return flag I then call the function constructor in my main, again like the tutorials: inner_region = Function('inner_region', inner_region_, extra_args=magnet) and define the region off the domain: magnet_region = domain.create_region('magnet0', 'vertices by inner_region','facet') (i think i also need to make this region a child of my Omega (all node) region) which fails miserably, due to the Function constructor seemingly not creating a Functions instance when called, but this is just from a bit of debugging (im no expert :D) Secondly, I can't seem to configure having two materials over differing regions in the mesh. I think this is partly due to the above but also very limited examples concerning this use case interactively. If someone could run me through this i'd be extremely grateful! BTW Im making meshes with pymesh, so cant use groups as a work around!

1 week, 5 days

ANN: SfePy 2021.2

by Robert Cimrman

I am pleased to announce the release of SfePy 2021.2. Description ----------- SfePy (simple finite elements in Python) is a software for solving systems of coupled partial differential equations by finite element methods. It is distributed under the new BSD license. Home page: https://sfepy.org Mailing list: https://mail.python.org/mm3/mailman3/lists/sfepy.python.org/ Git (source) repository, issue tracker: https://github.com/sfepy/sfepy Highlights of this release -------------------------- - new sensitivity analysis terms - positive FE basis based on Bernstein polynomials - smaller memory footprint of terms with constant material parameters For full release notes see [1]. Cheers, Robert Cimrman [1] http://docs.sfepy.org/doc/release_notes.html#id1 --- Contributors to this release in alphabetical order: Robert Cimrman Vladimir Lukes

2 weeks

Element stiffness matrix

by 979439736＠qq.com

I would like to ask what is the law of arrangement of the stiffness matrix of the element? I look forward to your reply, thank you!

2 weeks

Material coefficients with time derivatives of field variables

by Max Maximilian

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

2 weeks, 3 days

Implementation of IGA surface terms

by hugo.levy＠gmail.com

Hi Robert, I have been using Sfepy for the past two months now, and the more I use it the more I like it! Thank you for developing this great FEM (& more) software. I spent some time trying the various functionalities sfepy offers and recently discovered the IGA mode, which I managed to use interactively. After some testing, it turns out isogeometric analysis might be well suited for my applications, where geometry plays a crucial role. So far, my IGA simulations outperform FEM ones, which really sparked my interest for this type of analysis. Having no background in IGA, I read references [1] and [2] to grasp the mathematical concepts at stake, and reference [3] for the Bézier extraction method. That being said, user's guide mentions that surface terms in IGA are not implemented yet in Sfepy, and I was wondering why? For instance, I would like to compute integrals over (D-1)-dimensional spaces (where D denotes the problem spatial dimension), but as expected, running the line: >>> pb.evaluate('ev_surface_integrate.3.Gamma(u)', mode='eval') results with >>> ValueError: unknown dof connectivity type! (surface) This error is raised in sfepy\discrete\iga\fields.py function setup_extra_data. Of course, I could get around this missing functionality by sampling the scalar field of interest over the desired surface/curve manually (i.e. with field.evaluate_at) and use scipy.integrate. But that would mean giving up on the 'exact' geometry description, which I am not keen on. Is there a particular reason for not implementing those surface terms in IGA, like a technical difficulty that I am not yet aware of? Otherwise, could you give me insights into where to start if I were to implement my own quadrature rule for computing surface integrals in IGA mode? I haven't figured out in which module the evaluation of Bernstein basis functions at Gauss quadrature points occurs. I know juggling with NURBS, Bézier meshes, etc. doesn't seem particularly easy. Yet I'm willing to give it a try. Thank you very much for your help. Best regards, Hugo L. ------------------ [1] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, Elsevier, 2005, 194 (39-41), pp.4135-4195. 10.1016/j.cma.2004.10.008. hal-01513346 [2] Cimrman, Robert. (2014). Enhancing SfePy with Isogeometric Analysis. [3] Michael J. Borden, Michael A. Scott, John A. Evans, Thomas J. R. Hughes: Isogeometric finite element data structures based on Bezier extraction of NURBS, Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, Texas, March 2010.

1 month

AttributeError: 'float' object has no attribute 'ndim'

by Max Maximilian

Hey all, I'm struggled with an issue called AttributeError: 'float' object has no attribute 'ndim'. It appears when starting at the initial condition calculation. Is there any problems with my material setting? #codes: 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 = 1001 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.01 porosity = 0.3 # PTL porosity theta_PTL = nm.pi*30./180. # PTL contact angle KG = 4.5e-14 # Permeability, m2 R = 8.3145 # Ideal gas constant, J/mol/K T = 273.15+80 # Ambient temperature, K patm = 101325.0 # Ambient pressure s0 = 0.5 # Channel liquid saturation Mox = 0.032 # Molar mass density of oxygen, kg/mol # Boundary pressure 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-s0)-2.120*(1-s0)**2+1.263*(1-s0)**3) rho00 = (patm+pc0)/R/T*Mox example_name = "fluid_dynamics_1D" def post_process(out, pb, state, extend=False): # Gradient evalution rhow_grad = pb.evaluate('ev_grad.4.Omega(rho)', mode='el_avg', verbose=False) out['pressure_gradient'] = Struct(name='output_data', mode='cell', data=rhow_grad, 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. """ porosity = 0.3 # PTL porosity theta_PTL = nm.pi*30./180. # PTL contact angle KG = 4.5e-14 # Permeability, m2 muG = 1.881e-5 # Gas viscosity, Pa*s muL = 3.56e-4 # Liqiud viscosity, Pa*s R = 8.3145 # Ideal gas constant, J/mol/K T = 273.15+80 # Ambient temperature, K Mox = 0.032 # Molar mass density, kg/mol rho_H2O = 996 # Water density, kg/m3 m = 3 # Relative permeabiltiy coefficient Mw = 0.018 # Molar mass density of water, kg/mol i_app = 12000. # Applied current density, A/m2 F = 96485 # Faraday constant, C/mol gas_rate = -i_app/4/F*Mox liquid_rate = i_app/2/F*Mw dpcds = nm.zeros(coors.shape[0]) if mode == 'qp': # Field and gradient values in quadrature points coordinates given by integral i # - they are the same as in `coors` argument. liquid_saturation_values = problem.evaluate('ev_volume_integrate.i.Omega(s)', mode='qp', verbose=False) liquid_saturation_values.shape = (coors.shape[0], 1, 1) 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) grad_rho_values = problem.evaluate('ev_grad.i.Omega(rho)', mode='qp', verbose=False) grad_rho_values.shape = (coors.shape[0], 1, 1) rho_values = problem.evaluate('ev_volume_integrate.i.Omega(rho)', mode='qp', verbose=False) rho_values.shape = (coors.shape[0], 1, 1) # 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 # Stiffness matrix of gas density: rho*KG*(1-s)^3/muG*R*T/Mox stiff_dens = rho_values*KG*(1-liquid_saturation_values)**m/muG*R*T/Mox # 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) # Diversity matrix of liquid saturation: rho_H2O*K*s^m/muL*R*T/Mox*grad_rho div_sat = rho_H2O*KG*liquid_saturation_values**m/muL*R*T/Mox*grad_rho_values # Stiffness matrix of liquid saturation: rho_H2O*dpcds 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*dpcds return {'mass_dens' : mass_dens, 'source_dens' : source_dens, 'stiff_dens' : stiff_dens, 'mass_sat' : mass_sat, 'div_sat' : div_sat, 'stiff_sat' : stiff_sat, 'gas_flux' : gas_rate, 'liquid_flux' : liquid_rate} materials = { 'coef' : 'get_coef', } fields = { 'gas_density': ('real', 'scalar', 'Omega', 2), 'liquid_saturation': ('real', 'scalar', 'Omega', 2), } 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'), } regions = { 'Omega' : 'all', 'Gamma_Left' : ('vertices in (x < 1e-8)', 'facet'), 'Gamma_Right' : ('vertices in (x > 1.99999e-4)', 'facet'), } ebcs = { 'channel_dens' : ('Gamma_Right', {'rho.0' : rho00}), 'channel_satu' : ('Gamma_Right', {'s.0' : s0}), } 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_laplace.i.Omega(coef.stiff_dens, v1, rho) = dw_surface_integrate.i.Gamma_Left(coef.gas_flux, v1)""", 'liquid_saturation' : """dw_volume_dot.i.Omega(coef.mass_sat, v2, ds/dt) + dw_diffusion_r.i.Omega(coef.div_sat, v2) - dw_laplace.i.Omega(coef.stiff_sat, v2, s) = -dw_surface_integrate.i.Gamma_Left(coef.liquid_flux, v2)""", } functions = { 'get_coef' : (get_coef,), } ics = { 'ic' : ('Omega', {'rho.0' : rho00, 's.0' : s0}), } solvers = { 'ls' : ('ls.scipy_direct', { }), 'newton' : ('nls.newton', { 'eps_a' : 1e-8, 'eps_r' : 1e-2, 'problem' : 'nonlinear', }), 'tss' : ('ts.adaptive', { 't0' : t0, 't1' : t1, 'dt' : None, 'verbose' : 1, 'n_step' : 100, 'quasistatic' : True, }), } options = { 'ts' : 'tss', 'nls' : 'newton', 'ls' : 'ls', 'save_times' : 10, 'output_dir' : 'output/' + example_name, 'output_format' : "vtk", 'post_process_hook' : 'post_process', }

1 month, 4 weeks

Elastodynamic calculation in interactive mode

by fabien.tholin＠onera.fr

Hi everyone. I am using Sfepy for several weeks and i found it amazing. I have some trouble to use it interactively: I would like to run an elastodynamic calculation in interactive mode, but my code seems to give the steady state solution at each output and i really do not know what is missing. The configuration is a 2D plate subject to a pressure field on its top surface. I would like to see the dynamics of the deformation on short timescales. Probably it is a basic mistake but i am really stuck with it. Thank you very much for your help. Regards. Here is my code: #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Fri May 7 15:40:38 2021 2D elastodynamics test-case: A steel plate is subjet to a Gaussian pressure field centered on its top boundary """ from os.path import join as pjoin import numpy as nm from sfepy.discrete.fem import FEDomain, Field from sfepy.mesh.mesh_generators import gen_block_mesh from sfepy.discrete import (FieldVariable, Material, Integral, Equation, Equations, Problem) #from sfepy.mechanics.matcoefs import stiffness_from_lame import sfepy.mechanics.matcoefs as mc from sfepy.terms import Term from sfepy.discrete.conditions import Conditions, EssentialBC from sfepy.base.base import IndexedStruct from sfepy.solvers.ls import ScipyDirect from sfepy.solvers.nls import Newton from sfepy.solvers.ts_solvers import SimpleTimeSteppingSolver from sfepy.base.base import Struct outputs_folder = "out_sfepy" problem_name = "elastodynamics" output_folder = pjoin(outputs_folder, problem_name) output_format = "vtk" #================================= Mesh ====================================== #----------- generation of a 2D block mesh: Ndim=2 Ly = 2.0e-3 ; Lx = 10.0e-2 dims = [Lx, Ly] shape = [1000, 20] # Element size. HX = Lx / (shape[0] - 1) HY = Ly / (shape[1] - 1) # generation of the mesh mesh = gen_block_mesh(dims, shape, 0.5 * nm.array(dims), name='user_block', verbose=False) #=========================== Create a domain ================================= domain = FEDomain('domain', mesh) #=============== Define volume "Omega" and boundaries "Gamma" ================ min_x, max_x = domain.get_mesh_bounding_box()[:, 0] min_y, max_y = domain.get_mesh_bounding_box()[:, 1] eps = 1e-8 * (max_x - min_x) omega = domain.create_region('Omega', 'all') gamma1 = domain.create_region('Gamma1','vertices in x < %.10f' % (min_x + eps),'facet') gamma2 = domain.create_region('Gamma2','vertices in x > %.10f' % (max_x - eps),'facet') gamma3 = domain.create_region('Gamma3','vertices in y > %.10f' % (max_y - eps),'facet') gamma4 = domain.create_region('Gamma4','vertices in y < %.10f' % (min_y + eps),'facet') #================= Define the finite element approximation =================== # Using the field fu, we can define both the unknown variable and the test variable . field = Field.from_args('fu', nm.float64, 'vector', omega, approx_order=2) u = FieldVariable('u', 'unknown', field) du = FieldVariable('du', 'unknown', field) ddu = FieldVariable('ddu', 'unknown', field) v = FieldVariable('v', 'test', field, primary_var_name='u') dv = FieldVariable('dv', 'test', field, primary_var_name='du') ddv = FieldVariable('ddv', 'test', field, primary_var_name='ddu') #============================ Material parameters ============================ # from E, nu : E = 200.0e9 ; nu = 0.3 lam, mu = mc.lame_from_youngpoisson(E, nu, plane='strain') #============================ For elastodynamics ============================= rho = 7800.0 # Longitudinal and shear wave propagation speeds. cl = nm.sqrt((lam + 2.0 * mu) / rho) cs = nm.sqrt(mu / rho) #============================= Stiffness matrix ============================== m = Material('m', D=mc.stiffness_from_lame(dim=2, lam=lam, mu=mu),values={'rho':rho}) #================================= Time ====================================== # Time-stepping parameters. # Note: the Courant number C0 = dt * cl / H dt = HX / cl # C0 = 1 print('dt=', dt) timef = 1.5 * Lx / cl #======================== Pressure distribution ============================== def test_P(x, x0=0.0,R=1.0,P0=1.0): #--- champ de pression sinusoidal: #return P0*nm.sin(2.0*nm.pi*(x-x0)/R) #--- pic de pression return P0*nm.exp(-((x-x0)/R)**2.0) def tractionLoad( ts, coors, mode='qp', equations=None, term=None, problem=None): """ts = TimeStepper, coor = coordinates of field nodes in region.""" if mode == 'qp': val=nm.zeros((coors.shape[0], 3, 1),dtype=nm.float64) for i in range(coors.shape[0]): x=coors[i,0] #; y=coors[i,1] pressure=test_P(x,0.5*dims[0], dims[0], 1.0e7) val[i,:,0]=0.0, pressure, 0.0 # S11, S22, S12 return {'val' : val} p = Material(name='load',kind='time-dependent',function=tractionLoad) #=============================== Integration method ========================== integral = Integral('i', order=3) #========================== Definition of the "terms" ======================== t1 = Term.new('dw_lin_elastic(m.D, v, u)',integral, omega, m=m, v=v, u=u) t2 = Term.new('dw_volume_dot(m.rho, ddv, ddu)',integral, omega, m=m, ddv=ddv, ddu=ddu) #t3 = Term.new('dw_zero(dv, du)',integral, omega, m=m, dv=dv, du=du) t4 = Term.new('dw_surface_ltr(load.val, v )',integral, gamma3, load=p,v=v) eq = Equation('balance_of_forces in time',t1+t2+t4) eqs = Equations([eq]) #============================ boundary conditions ============================ fix_u1 = EssentialBC('fix_u', gamma1, {'u.all' : 0.0}) fix_u2 = EssentialBC('fix_u', gamma2, {'u.all' : 0.0}) #=========================== Problem instance ================================ pb = Problem(problem_name,equations=eqs) pb.set_bcs(ebcs=Conditions([fix_u1, fix_u2])) pb.setup_output(output_dir=output_folder, output_format=output_format) pb.update_materials() #=============================== Solveur ===================================== ls = ScipyDirect({'use_presolve' : True}) nls_status = IndexedStruct() nls = Newton({'i_max': 1, 'eps_a': 1e-6, 'eps_r' : 1e-6,}, lin_solver=ls, status=nls_status) tss = SimpleTimeSteppingSolver({'t0' : 0.0, 't1' : 1.0e-9, 'n_step' : 10}, nls=nls, context=pb, verbose=True) #================================ SOLVE ====================================== pb.set_solver(tss) status = IndexedStruct() vec = pb.solve(status=status) print('Nonlinear solver status:\n', nls_status) print('Stationary solver status:\n', status) pb.save_state('linear_elasticity.vtk', vec)

1 month, 4 weeks

1D problem within the interactive way

by Max Maximilian

2 months, 1 week

Material Parameters

by 979439736＠qq.com

I want different units to set different elastic modulus. What should I do？Thanks！

2 months, 1 week

Jump to page:

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

SfePy