Hi Fabien,
On 28. 05. 21 21:05, fabien.tholin@onera.fr wrote:
Hi everyone.
I am using Sfepy for several weeks and i found it amazing.
Thanks!
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.
The main issue is that you need to use one of the elastodynamic solvers (e.g. ts.bathe). More comments below in the code listing.
Thank you very much for your help.
Regards.
Here is my code:
<snip>
#=============================== Integration method ==========================
integral = Integral('i', order=3)
You may need quadrature order=4, but it seems to be working with 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)
The elastodynamic solvers require including the velocity-dependent term even if it is zero, so that the matrices are correctly allocated:
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+t3+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})
Also du, ddu need to be fixed:
fix_u1 = EssentialBC('fix_u', gamma1, {'u.all' : 0.0, 'du.all' : 0.0, 'ddu.all' : 0.0}) fix_u2 = EssentialBC('fix_u', gamma2, {'u.all' : 0.0, 'du.all' : 0.0, 'ddu.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)
Replace tss for example with:
from sfepy.solvers import Solver tss = Solver.any_from_conf( Struct(**{ 'kind' : 'ts.bathe', 't0' : 0.0, 't1' : 1.0e-6, 'dt' : dt, 'is_linear' : True, 'verbose' : 1, }), nls=nls, context=pb, )
or other tss listed in examples/linear_elasticity/elastodynamic.py.
r.