I want to find the reaction force at the end of the total dynamic displacement given. There is no dynamic force given to the system. Only the displacement is provided at the boundary of a plate in some steps. Displacement is not fixed at each step which makes it a dynamic system. As I am considering the dynamic case so acceleration should also be present in the weak form. But, I am not able to consider the effect of that acceleration of the previous step, when I am trying to solve displacement and reaction at the next step. Please help…
This is the code i have made.
from dolfin import *
import fenics as fe
from fenics import *
import matplotlib.pyplot as plt
from ufl import nabla_grad
from ufl import nabla_div
import numpy as np
import pygmsh
import gmsh
import meshio
import pandas as pd
import os
import sympy as sp
import math
loadsteps = list(np.arange(0,61,10))
loadsteps = [10]
No_of_timestep = 50
d_t = 1/No_of_timestep
dt = Constant(d_t)
time = np.array(np.arange(dt,1+2*dt,dt))
u_x_right = np.zeros(No_of_timestep+1)
u_y_top = np.zeros(No_of_timestep+1)
Generalized-alpha method parameters
alpha_m = Constant(0.2)
alpha_f = Constant(0.4)
gamma = Constant(0.5+alpha_f-alpha_m)
beta = Constant((gamma+0.5)**2/4.)
del_t = np.array(np.arange(0,1+d_t,d_t))
delta = 0.1
for l in range(0,No_of_timestep+1):
if (l==0):
u_x_right[l] = 0
u_y_top[l] = 0
else:
u_x_right[l] = 10 + math.sin(21math.pidel_t[l])+math.cos(11math.pidel_t[l])+math.sin(1math.pidel_t[l])
u_y_top[l] = u_x_right[l]
u_x_right = (0.05/sum(u_x_right))*u_x_right
u_y_top = (0.1/sum(u_y_top))*u_y_top
u_x_right_bc = u_x_right.cumsum()
u_y_top_bc = u_y_top.cumsum()
np.delete(u_x_right, 0)
np.delete(u_y_top, 0)
Rxlb = np.zeros(No_of_timestep+1)
Rxrb = np.zeros(No_of_timestep+1)
Rybb = np.zeros(No_of_timestep+1)
Rytb = np.zeros(No_of_timestep+1)
Displacement in x direction in steps
epsilon_11 = np.zeros(len(u_x_right))
sigma_11 = np.zeros(len(epsilon_11))
Sig11_anly = np.zeros(len(epsilon_11))
Displacement in y direction in steps
epsilon_22 = np.zeros(len(u_y_top))
sigma_22 = np.zeros(len(epsilon_22))
Sig22_anly = np.zeros(len(epsilon_22))
Mass density
rho = Constant(1.0)
Mass form
def m(a, du):
return rho*inner(a, du)*dx
Update formula for acceleration
def update_a(u, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
beta_ = beta
else:
dt_ = float(dt)
beta_ = float(beta)
return (u-u_old-dt_v_old)/beta_/dt_**2 - (1-2beta_)/2/beta_*a_old
Update formula for velocity
def update_v(a, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
gamma_ = gamma
else:
dt_ = float(dt)
gamma_ = float(gamma)
return v_old + dt_*((1-gamma_)*a_old + gamma_a)
def avg(x_old, x_new, alpha):
return alphax_old + (1-alpha)*x_new
for j in range(0,len(loadsteps)):
print("loadstep = ",loadsteps[j])
## Dimensions
L_x = 1.0
L_y = 1.0
## Meshing
gmsh.initialize()
with pygmsh.occ.Geometry() as geom:
geom.characteristic_length_min = 0.0
geom.characteristic_length_max = 0.03
rectangle = geom.add_rectangle([0, 0, 0], L_x, L_y)
mesh = geom.generate_mesh()
triangle_cells = None
for cell in mesh.cells:
if cell.type == "triangle":
# Collect the individual meshes
if triangle_cells is None:
triangle_cells = cell.data
else:
triangle_cells = np.vstack((triangle_cells, cell.data))
triangle_data = None
triangle_mesh = meshio.Mesh(points=mesh.points, cells={"triangle": triangle_cells})
nodes1 = mesh.points[:,0:2]
cells = triangle_cells
nodes = nodes1
for t in range(1,No_of_timestep+1):
mesh = Mesh()
editor = MeshEditor()
editor.open(mesh,"triangle", 2, 2)
editor.init_vertices(len(nodes))
editor.init_cells(len(cells))
[editor.add_vertex(i,n) for i,n in enumerate(nodes)]
[editor.add_cell(i,n) for i,n in enumerate(cells)]
editor.close()
V = VectorFunctionSpace(mesh, "Lagrange", 1)
Du = TrialFunction(V)
du = TestFunction(V)
u = Function(V, name="displacement")
a = Function(V)
# if (t==1):
v = Function(V)
u_t_minus = Function(V)
a_t_minus = Function(V)
v_t_minus = Function(V)
I = Identity(len(u))
## Defining Boundaries
def bottom(x, on_boundary):
return (on_boundary and fe.near(x[1], 0.0))
def top(x, on_boundary):
return (on_boundary and fe.near(x[1], L_y+u_y_top_bc[t-1], 1.00e-5))
def left(x, on_boundary):
return (on_boundary and fe.near(x[0], 0.0))
def right(x, on_boundary):
return (on_boundary and fe.near(x[0], L_x+u_x_right_bc[t-1], 1.00e-5))
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
AutoSubDomain(left).mark(boundary_subdomains, 1)
AutoSubDomain(right).mark(boundary_subdomains, 2)
AutoSubDomain(bottom).mark(boundary_subdomains, 3)
AutoSubDomain(top).mark(boundary_subdomains, 4)
dss = ds(subdomain_data=boundary_subdomains)
## Dirichlet Boundary conditions
## For simply support ends
## left end restricted to move in x direction but free to move in y direction
## right end is given u_x_right displacement in x direction and free to move in y direction
bc_l = DirichletBC(V.sub(0), 0.0, left)
left_bdry_dofs = np.fromiter(bc_l.get_boundary_values().keys(), dtype = int)
bc_r = DirichletBC(V.sub(0), u_x_right[t], right)
right_bdry_dofs = np.fromiter(bc_r.get_boundary_values().keys(), dtype = int)
bc_b = DirichletBC(V.sub(1), 0.0, bottom)
bottom_bdry_dofs = np.fromiter(bc_b.get_boundary_values().keys(), dtype = int)
bc_t = DirichletBC(V.sub(1), u_y_top[t], top)
top_bdry_dofs = np.fromiter(bc_t.get_boundary_values().keys(), dtype = int)
bcs = [bc_l,bc_r,bc_b,bc_t]
# Kinematics
B = Constant((0, 0))
T = Constant((0, 0))
F = I + grad(u) # Deformation gradient
J = det(F)
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
I1 = tr(C) + 1
I2 = 0.5*((tr(C)+1)**2 - (tr(C*C)+1))
I3 = J**2
I1_bar = J**(-2/3)*I1
I2_bar = J**(-4/3)*I2
# Stored strain energy density
W = 0.5*(I1_bar - 3) + 1.5*((J-1)**2)
a = update_a(Du,u_t_minus, v_t_minus, a_t_minus, ufl = True)
f_int = m(avg(a_t_minus, a, alpha_m), du) + W*dx
f_ext = dot(B, u)*dx + dot(T, u)*dss
R = f_int - f_ext
F = derivative(R, u, du)
Jacobian = derivative(F, u, Du)
solve(F == 0, u, bcs, J=Jacobian)
u_x = np.zeros(len(nodes))
u_y = np.zeros(len(nodes))
for n in range(0,len(nodes)):
u_x[n] = u(nodes[n,0],nodes[n,1])[0]
u_y[n] = u(nodes[n,0],nodes[n,1])[1]
f_ext_unknown = assemble(F)
## Sorting x and y degrees of freedom
x_dofs = V.sub(0).dofmap().dofs()
y_dofs = V.sub(1).dofmap().dofs()
## Sorting the reactions in x and y direction residuals of boundaries
## x component of residuals of nodes present in at left boundary
Rxlb_step = [0 for element in range(len(left_bdry_dofs))]
k = 0
for i in left_bdry_dofs:
Rxlb_step[k] = f_ext_unknown[i]
k +=1
## x component of residuals of nodes present in right boundary
Rxrb_step = [0 for element in range(len(right_bdry_dofs))]
k = 0
for i in right_bdry_dofs:
Rxrb_step[k] = f_ext_unknown[i]
k +=1
## y component of residuals of nodes present in at bottom boundary
Rybb_step = [0 for element in range(len(bottom_bdry_dofs))]
k = 0
for i in bottom_bdry_dofs:
Rybb_step[k] = f_ext_unknown[i]
k +=1
## y component of residuals of nodes present in top boundary
Rytb_step = [0 for element in range(len(top_bdry_dofs))]
k = 0
for i in top_bdry_dofs:
Rytb_step[k] = f_ext_unknown[i]
k +=1
Rxlb[t] = sum(Rxlb_step)
Rxrb[t] = sum(Rxrb_step)
Rybb[t] = sum(Rybb_step)
Rytb[t] = sum(Rytb_step)
reactions = [sum(Rxlb),sum(Rxrb),sum(Rybb),sum(Rytb)]
# Get vectors (references)
u_vec, u_t_minus_vec = u.vector(), u_t_minus.vector()
v_t_minus_vec, a_t_minus_vec = v_t_minus.vector(), a_t_minus.vector()
a.vector()[:] = update_a(u_vec,u_t_minus_vec,v_t_minus_vec,a_t_minus_vec,ufl = True)
a_vec = a.vector()
v.vector()[:] = update_v(a_vec,v_t_minus_vec,a_t_minus_vec, ufl = True)
u_t_minus.vector()[:] = u.vector()
a_t_minus.vector()[:] = a.vector()
v_t_minus.vector()[:] = v.vector()
step_displacements = np.vstack((u_x, u_y)).T
nodes = nodes + step_displacements
## Plotting the final state
fe.plot(u, mode="displacement")
plt.show()
total_displacement = nodes - nodes1