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