Hello Everyone,
I am trying to include an extra energy term in the standard elastic strain energy term. The extra term depends on a function of strain and activates only after a certain condition. Addition of this term is not showing its effect on the results. According to my understanding: the extra term is not directly depends on the primary variable and vanishes during differentiation of strong form. I have attached the minimal code here. In the code the extra energy term is residual_energy
from fenics import *
import numpy as np
from mshr import *
import os, shutil
E = 24252.
nu = 0.2
mu = Expression('E_1/(2.0*(1.0 + nu_1))', E_1 = E, nu_1 = nu,degree = 0)
lmbda = Expression('(E_1 * nu_1)/((1.0 - 2.0*nu_1)*(1.0 + nu_1))', E_1 = E, nu_1 = nu,degree = 0)
G_strain = 3.756
H_strain = 1.1
domain = Rectangle(Point(0.0, 0.0), Point(100.0, 100.0))
mesh = generate_mesh(domain, 32)
ndim = mesh.topology().dim()
thickness = 50.
V_u = VectorFunctionSpace(mesh, 'CG', 1)
V_plot = FunctionSpace(mesh, 'DG', 0)
u, du, v = Function(V_u), TrialFunction(V_u), TestFunction(V_u)
zero_1 = Constant(0.0)
eps_p_0 = project(zero_1,V_plot)
eps_p1 = Expression('eps_p1',eps_p1 = eps_p_0,degree = 0)
u.rename('displacement','displacement')
tol = 1E-5
tol_v = 1e-20
def boundary_D_l(x, on_boundary):
return x[0] <= 0.0 + tol and x[1] <= 0.0 + tol
def boundary_D_r(x, on_boundary):
return x[0] <= 0.0 + tol
def boundary_D_t(x, on_boundary):
return x[1] >= 100.0 - tol
bcl = DirichletBC(V_u, Constant((0.0 ,0.0)), boundary_D_l, method='pointwise')
bcr = DirichletBC(V_u.sub(1), Constant(0.0), boundary_D_r, method='pointwise')
def eps(u,eps_p1):
mat = as_matrix([[1., 0.], [0., 0.]])
return sym(grad(u)) - eps_p1*mat
def eps_eig(eps_p1):
mat = as_matrix([[1., 0.], [0., 0.]])
return eps_p1*mat
def sigma(u):
return 2.0*mu*(eps(u,eps_p1)) + lmbda*tr(eps(u,eps_p1))*Identity(ndim)
num_steps = 5
total_disp = 0.5
disp_steps = total_disp/(num_steps)
elastic_energy_1 = 0.5*lmbda*(tr(eps(u,eps_p1)))**2 + mu*tr(eps(u,eps_p1)*eps(u,eps_p1))
residual_energy = 0.5*lmbda*(tr(eps_eig(eps_p1)))**2 + mu*tr(eps_eig(eps_p1)*eps_eig(eps_p1))
total_energy = (elastic_energy_1 + residual_energy)*thickness*dx
E_u = derivative(total_energy,u,v)
Jd = derivative(E_u, u, du)
u_R = Expression(('disp_app'),disp_app = disp_steps, degree=0)
bct = DirichletBC(V_u.sub(1), u_R, boundary_D_t, method='pointwise')
bc_disp = [bcl, bcr, bct]
problem_u = NonlinearVariationalProblem(E_u, u, bc_disp, Jd)
solver_u = NonlinearVariationalSolver(problem_u)
def alternate_minimization(u,n):
solver_u.solve()
if n >= 2:
sn_val = eps(u,eps_p1)
sn_c = (2/3)*sqrt(0.5*((sn_val[0,0]-sn_val[1,1])**2+(sn_val[1,1])**2+(sn_val[0,0])**2)+(3*(sn_val[0,1]**2)))
eps_u = project(sn_c,V_plot)
eps_p = project(G_strain*((eps_u)**H_strain),V_plot)
eps_p1.eps_p1 = eps_p
savedir = "Results/"
if os.path.isdir(savedir):
shutil.rmtree(savedir)
file_u = File(savedir+"/u.pvd")
forces = np.zeros((num_steps+1, 2))
def postprocessing():
b_int = assemble(E_u)
bcl.apply(b_int)
bcr.apply(b_int)
fint=b_int.copy()
forc = -fint[int(list(bct.get_boundary_values().keys())[1])]
forces[n+1] = np.array([u(100.,100.)[1],forc])
file_u << (u,n)
np.savetxt(savedir+'/forces.txt', forces)
for n in range(num_steps):
u_R.disp_app = disp_steps*(n+1)
alternate_minimization(u,n)
postprocessing()
How I can add the effect of this extra term in the solution. Please help to find the solution of this problem.
Thank you,
Manish