Updated lagrangian elastic problem

variational formulation
1

Hi,
I’m working on an updated lagrangian formulation of hyperelastic 2D problem with a phase field damage variable.

from dolfin import *
from mshr import*

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

snes_solver_parameters = {"nonlinear_solver": "snes",
                        "symmetric": True,
                        "snes_solver": {"linear_solver": "lu",
                                        "preconditioner": "petsc_amg",
                                        "maximum_iterations": 200,
                                        "report": True,
                                        "line_search": 'l2',
                                        "error_on_nonconvergence": True,
                                        "relative_tolerance": 1e-7, 
                                        "absolute_tolerance": 1e-7}}

#Geometry
n_mesh=100
domain=Rectangle(Point(0,0),Point(1,.1))
mesh0=generate_mesh(domain, n_mesh)
mesht=generate_mesh(domain, n_mesh)
type_mesh=''
tollb=1E-13
def boundary_left(x): return abs(x[0]) < tollb
def boundary_right(x):return abs(x[0]-1) < tollb
def boundary_bottom(x):return abs(x[1]) <tollb
#functional space
   #mesh0
Space1D1_0=FunctionSpace(mesh0,'P',1)
Space2D2_0=VectorFunctionSpace(mesh0,'P',1,2)
SpaceMatrixnonsym_0=TensorFunctionSpace(mesh0,'DG',0,(2,2))
   #mesht
Space1D1_t=FunctionSpace(mesht,'P',1)
Space2D2_t=VectorFunctionSpace(mesht,'P',1,2)
SpaceMatrixnonsym_t=TensorFunctionSpace(mesht,'DG',0,(2,2))

d_t=Function(Space1D1_t)
u_t=Function(Space2D2_t)
u_Test_t=TestFunction(Space2D2_t)
u_Trial_t=TrialFunction(Space2D2_t)

B_t=project(Identity(2),SpaceMatrixnonsym_t)

d_0=Function(Space1D1_0)
dpre_0=Function(Space1D1_0)
d_step_pre_0=Function(Space1D1_0)
Test_d_0=TestFunction(Space1D1_0)
Trial_d_0=TrialFunction(Space1D1_0)

u_0=Function(Space2D2_0)
u_global_0=Function(Space2D2_0)
delta_u_0=Function(Space2D2_0)

#parameters

#elastic

mu=[2.4,.008]
alpha=[.5,2.8]
m_1=[2.4,.008]
alpha_1=[.5,2.8]

#damage
lam_crack=3

Du_damage=.1

I1_crack=(lam_crack**2)+(2*(lam_crack**(-1)))
psi_crack=sum([(mu[a]*((3**(1-alpha[a])))/(2*alpha[a]))*((I1_crack**alpha[a]) -(3**alpha[a])) for a in range(len(mu))])
G_frac=4*psi_crack*Du_damage/3

b_int=Du_damage/4
c_0_dam=8/3


#evolution
dotu=1
dt=.01
tfin=1.2*(lam_crack-1)
t=0

#functions

def F(u): return grad(u) + Identity(2)
def C(u):return F(u).T*F(u)
def B(u): return F(u)*(F(u).T)

def I_1_global(u): return tr(C(u))+(det(C(u))**(-1))

def degr_func_inf(d): return  (1-d)**2
def degr_func_1(d): return   (1-d)**2
def alpha_dam(d): return d
def gamma_dam(d): return (1/(c_0_dam*b_int))*((alpha_dam(d))+((b_int**2)*dot(grad(d),grad(d))))


def elastic_energy_infinity_0(u):
    return sum([(mu[a]*((3**(1-alpha[a])))/(2*alpha[a]))*((I_1_global(u)**alpha[a]) -(3**alpha[a])) for a in range(len(mu))])

def energy_functional_0(u,d): return (degr_func_inf(d)*elastic_energy_infinity_0(u))+(G_frac*(gamma_dam(d)))

def elastic_energy_infinity_t(Bt,ut): 
    return sum([(mu[a]*((3**(1-alpha[a])))/(2*alpha[a]))*(((inner(Bt,C(ut))+((det(Bt)*det(C(ut)))**(-1)))**alpha[a])-(3**alpha[a])) for a in range(len(mu))])

def energy_functional_t(Bt,ut,d):return degr_func_inf(d)*elastic_energy_infinity_t(Bt,ut)

def stress_func(u,d):
    req=  sum([(mu[a]*( 3**(1-alpha[a])))*((I_1_global(u)**(alpha[a]-1))) for a in range(len(mu))])
    return degr_func_inf(d)*req*(F(u)-((1/det(C(u)))*inv((F(u)).T)))

#problem for damage

d_min=interpolate(Constant(0),Space1D1_0)
d_max=interpolate(Constant(1-(1e-3)),Space1D1_0)
#boundary conditions Dirichlet (5 lines), if you want use the Newmann b.c., comment the next 5 lines
bcd_dleft=DirichletBC(Space1D1_0,0,boundary_left)
bc_dright=DirichletBC(Space1D1_0,0,boundary_right)
bcd=[bcd_dleft,bc_dright]
for bc in bcd:
    bc.apply(d_max.vector())
class DamageProblem(OptimisationProblem):
    def f(self,x):
        d_0.vector()[:]=x
        return assemble(energy_functional_0(u_global_0,d_0)*dx(mesh0))
    def F(self,b,x):
        d_0.vector()[:]=x
        assemble(derivative(energy_functional_0(u_global_0,d_0)*dx(mesh0),d_0,Test_d_0),b)
    def J(self,A,x):
        d_0.vector()[:]=x
        assemble(derivative(derivative(energy_functional_0(u_global_0,d_0)*dx(mesh0),d_0,Test_d_0),d_0,Trial_d_0),A)
solver_d=PETScTAOSolver()
solver_d.parameters.update({'maximum_iterations':400,"gradient_relative_tol" : 1.0e-4,'report':True})

tol=1e-7
maxite=20
while t<tfin:
    err=1
    ite=0
    x, y = mesht.coordinates()[:, 0], mesht.coordinates()[:, 1]
    lpxt=max(x)
    while err>tol and ite<maxite:
        d_t.vector()[:]=d_0.vector()[:]
        x, y = mesht.coordinates()[:, 0], mesht.coordinates()[:, 1]
        lpxt=max(x)
        def boundary_left_t(x): return abs(x[0]) < tollb
        def boundary_right_t(x):return abs(x[0]-lpxt) < tollb
        def boundary_bottom_t(x):return abs(x[1]) <tollb
        bc_left=DirichletBC(Space2D2_t, Constant((0,0)), boundary_left_t)
        #bc_right=DirichletBC(Space2D2_t,Constant((dimp,0)),boundary_right_t)
        if ite==0:
            bc_right=DirichletBC(Space2D2_t,Constant((dt*dotu,0)),boundary_right_t)
        else:
            bc_right=DirichletBC(Space2D2_t,Constant((0,0)),boundary_right_t)
        bc=[bc_left,bc_right]
        energy=energy_functional_t(B_t,u_t,d_t)*dx(mesht)
        var1=derivative(energy,u_t,u_Test_t)
        var2=derivative(var1,u_t,u_Trial_t)
        problem = NonlinearVariationalProblem(var1,u_t,bc,var2)
        solver = NonlinearVariationalSolver(problem)
        solver.parameters.update(snes_solver_parameters)
        solver.solve()
        delta_u_0.vector()[:]=u_t.vector()[:]
        u_supp=project(u_0+delta_u_0,Space2D2_0)
        u_global_0.assign(u_supp)#(u_0+delta_u_0)
        solver_d.solve(DamageProblem(),d_0.vector(),d_min.vector(),d_max.vector())
        deltad=project(d_0-dpre_0,Space1D1_0)
        #update mesh
        u_0=project(u_global_0,Space2D2_0) 
        for point_v in range(len(mesh0.coordinates())):
            mesht.coordinates()[point_v]=mesh0.coordinates()[point_v]+u_0(Point(mesh0.coordinates()[point_v]))
        Space1D1_t=FunctionSpace(mesht,'P',1)
        Space2D2_t=VectorFunctionSpace(mesht,'P',1,2)
        SpaceMatrixnonsym_t=TensorFunctionSpace(mesht,'DG',0,(2,2))
        d_t=Function(Space1D1_t)
        u_t=Function(Space2D2_t)
        u_Test_t=TestFunction(Space2D2_t)
        u_Trial_t=TrialFunction(Space2D2_t)   
        B_t_support=project(Identity(2),SpaceMatrixnonsym_t)
        B_0=project(B(u_0),SpaceMatrixnonsym_0)
        B_t_support.vector()[:]=B_0.vector()[:]
        B_t=project(B_t_support,SpaceMatrixnonsym_t)
        err=max(deltad.vector())
        dpre_0.assign(d_0)
        ite=ite+1

    #update parameters
    u_0=project(u_global_0,Space2D2_0) 
    d_min.assign(d_0)
    t=t+dt

The problem is that the evaluation of gradient of u_0 (and then of B_0 and B_t) is not precisely and, thus the yield condition is reached in a wrong way. How I can solve this isuue?

Thanks for attention and excuseme for the long topic

1 Like
2

Please use markdown syntax to encapsulate code, ie.

```python
def f(x):
    return x[0]
```

and make sure that your code can be executed with copy-pasting it into an editor.

Also your code is rather long and complicated. To maximize the likelihood of getting a response, I would suggest that you try to make a small as possible example to reproduce the issue.

3

thanks for the answer, it’s ok now? I’m trying to reduce it more