Thank you very much for your help ! It is true that i completely forgot to set the boundary conditions for time derivatives.
Now it seems to work perfectly. Now this script compute the dynamics of the plate due to a pressure distribution varying in space and time. I put here my last script if someone is interested:

!/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 #from examples.linear_elasticity.its2D_2 import stress_strain

output_folder = "out_sfepy" problem_name = "elastodynamics" output_format = "vtk" verb=False

#================================== Time =====================================

Noutput=100 dtoutput=1.0e-6 dtmin=1.0e-7 CFL=False

#================================= Mesh ======================================

#----------- generation of a 2D block mesh: Ndim=2 Ly = 2.0e-3 ; Lx=0.04 dims = [Lx, Ly] #shape = [400, 20] shape = [200, 10] # 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 = 20.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 dtcfl = HX / cl # C0 = 1 print('dtcfl=', dtcfl)

if(CFL): dt=nm.min([dtcfl,dtmin,dtoutput]) else: dt=nm.min([dtmin,dtoutput])

#======================== 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 biexp(t,wave='D',fac=1.0,norm=False):

if(wave == 'D'):
A=109405.0 ; alpha=22708.0 ; beta=1294530.0 ; Imax=100000.0 elif(wave == 'A'): A=218810.0 ; alpha=11354.0 ; beta=647265.0 ; Imax=200000.0 elif(wave == 'G'): A=-1044746.2664139411842 ; alpha=81169.502440601603134 ; beta=1.0/(0.08*200.0e-6) ; Imax=100000.0

if(norm): return fac*A*( nm.exp(-alpha*t) - nm.exp(-beta*t) )/Imax else: return fac*A*( nm.exp(-alpha*t) - nm.exp(-beta*t) )

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.0, 0.01, 1.0e8) val[i,:,0]=0.0, pressure, 0.0 # S11, S22, S12
t=ts.time val=val*biexp(t,'D',fac=1.0, norm=True) return {'val' : val}

p = Material(name='load',kind='time-dependent',function=tractionLoad)

#=============================== Integration method ==========================

integral = Integral('i', order=4)

#========================== 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)

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}) #fix_u1 = EssentialBC('fix_u', gamma1, {'u.all' : 0.0, 'du.all' : 0.0, 'ddu.all' : 0.0})

fix_u1 = EssentialBC('fix_u', gamma1, {'u.0' : 0.0, 'du.0' : 0.0, 'ddu.0' : 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])) #================================= Outputs ===================================

#=============================== 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 are listed in examples/linear_elasticity/elastodynamic.py. from sfepy.solvers import Solver

# tss = SimpleTimeSteppingSolver({'t0' : 0.0, 't1' : 1.0e-9, 'n_step' : 10}, nls=nls, context=pb, verbose=True)

#================================ SOLVE ====================================== t=0.0

for i in range(Noutput):

print('')
print('====================  iteration:', i , '===========================')  
print('')
print('    t=', t, ' t+dt=', t+dtoutput)

tss = Solver.any_from_conf(
    Struct(**{
        'kind' : 'ts.bathe',
        't0' : t, 't1' : t+dtoutput, 
        'n_step' : None, # precedence over dt
        'dt' : dt,
        'is_linear' : True,
        'verbose' : 0,
    }), nls=nls, context=pb,
)

pb.set_solver(tss)
status = IndexedStruct()    
#pb.update_materials()

if(i == 0):
    vec = pb.solve(state0=None,status=status,save_results=False)
else:
    vec = pb.solve(state0=vec,status=status,save_results=False)

if(verb):
    print('Nonlinear solver status:\n', nls_status)
    print('Stationary solver status:\n', status)

# ---Postprocess the solution. def vonmises_from_stress(stress): N=len(stress) vM=nm.zeros((N,1,1,1)) for i in range(N): #--- 2D: s11=stress[i,0,0,0] ; s22=stress[i,0,1,0] ; s12=stress[i,0,2,0] s33=0.0 ; s13=0.0 ; s23=0.0
#--- 3D: #s11=stress[i,0,0,0] ; s22=stress[i,0,1,0] ; s33=stress[i,0,2,0] #s12=stress[i,0,3,0] ; s13=stress[i,0,4,0] ; s23=stress[i,0,5,0]
vM[i,0,0,0]=nm.sqrt(0.5*((s11-s22)**2+(s22-s33)**2+(s33-s11)**2)+3.0*(s12**2+s23**2+s13**2))

    return vM

def stress_strain(out, pb, state, extend=False):
    """
    Calculate and output strain and stress for given displacements.
    """
    from sfepy.base.base import Struct

    ev = pb.evaluate
    strain = ev('ev_cauchy_strain.2.Omega(u)', mode='el_avg')
    stress = ev('ev_cauchy_stress.2.Omega(m.D, u)', mode='el_avg',
                copy_materials=False)
    
    
    vonmises=vonmises_from_stress(stress)

    out['cauchy_strain'] = Struct(name='output_data', mode='cell',
                                  data=strain, dofs=None)
    out['cauchy_stress'] = Struct(name='output_data', mode='cell',
                                  data=stress, dofs=None)
    
    out['vonMises']      = Struct(name='output_data', mode='cell',
                                  data=vonmises, dofs=None)

    return out

out = vec.create_output_dict()
out = stress_strain(out, pb, vec, extend=True)
#pb.setup_output(output_dir=output_folder, output_format=output_format)
pb.save_state(output_folder+'/'+problem_name+'_'+str(i)+'.'+output_format, out=out)

t=t+dtoutput