Hello,
I am trying to do an impact analysis on a plate for a nonlinear material based on Elastodynamics_Tutorial and Transient_heat_equation_using_mgis_fenics. I am using generalized alpha method based on the elastodynamics tutorial. The simulation runs without any issues for t>1s for cyclic loading cases, so the implementation seems to be right.
Since the simulation time for impact analysis is roughly 1e-3s, the size of the time increment is of the order 1e-4s~1e-5s. The following image shows the equation to compute the acceleration
Since my time increment is very small, the acceleration tends to be a large value because of the dt^2 term in the denominator.
Ideally, an explicit time integration should be used for dynamic analyses but I am using an implicit method since I use mgis_fenics which uses a NewtonSolver.
Is there a way to avoid the large accelerations using the same implicit scheme?
Please find the parts related to the dynamic implementation below:
def m(self, u, u_):
return self.rho*inner(u, u_)*self.problem.dx
def k(self, u_):
return (sum([inner(g.variation(u_), f.function) for (f, g) in zip(self.problem.fluxes.values(), self.problem.gradients.values())])*self.problem.dx)
def c(self, u, u_):
return self.eta_m*self.m(u, u_) + self.eta_k*self.k(u_)
def Wext(self, u_):
return self.loading*dot(self.n, u_)*self.ds(self.right_edge)
def update_a(self, u, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = self.dt
beta_ = self.beta
else:
dt_ = float(self.dt)
beta_ = float(self.beta)
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old
def update_v(self, a, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = self.dt
gamma_ = self.gamma
else:
dt_ = float(self.dt)
gamma_ = float(self.gamma)
return v_old + dt_*((1-gamma_)*a_old + gamma_*a)
def update_fields(self, u, u_old, v_old, a_old):
u_vec, u0_vec = u.vector(), u_old.vector()
v0_vec, a0_vec = v_old.vector(), a_old.vector()
a_vec = self.update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False)
v_vec = self.update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False)
print(v_vec.get_local(), a_vec.get_local())
v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec
u_old.vector()[:] = u.vector()
def avg(self, x_old, x_new, alpha):
return alpha*x_old + (1-alpha)*x_new
def compute_residual(self):
a_new = self.update_a(self.du, self.u_old, self.v_old, self.a_old, ufl=True)
v_new = self.update_v(a_new, self.u_old, self.v_old, self.a_old, ufl=True)
print("Residual", v_new, "\n", a_new)
residual = self.m(self.avg(self.a_old, a_new, self.alpha_m), self.u_) + self.c(self.avg(self.v_old, v_new, self.alpha_f), self.u_) \
+ self.k(self.u_)
if not self.disp_load:
residual -= self.Wext(self.u_)
return residual