# Phase field modeling

Hi everyone. I am doing a phase field modeling for fracture in julia using gridap. I have done two different method:1- is the phase field modeling by applying displacement on a boundary incrementally which works. 2- i am trying to do the same but not defining displacement on boundary. I have defined force implemented on boundary incrementally which is not working. I couldn’t understand why it is not working as every thing is the same. could you please tell me if you have any recommendation. thank you.

``````function   stepPhaseField(fem_params, uh_in,ΨPlusPrev_in)
A_PF(s,φ,uh_in,ΨPlusPrev_in) =Gc*ls*∇(φ)⋅∇(s)+ (2*ΨPlusPrev_in*s*φ) + (Gc/ls)*s*φ
a_PF(s,φ) = ∫(A_PF(s,φ,uh_in,ΨPlusPrev_in))fem_params.dΩ
b_PF(φ) =∫((Gc/ls)*φ)fem_params.dΩ
op_PF = AffineFEOperator(a_PF,b_PF ,U_PF ,V0_PF)
sh_out = solve(op_PF)
return  get_free_dof_values(sh_out)
end
function  stepDisp(fem_params,pth, uh_in,sh_in, vApp)
g = VectorValue(0, -vApp)
A_Disp(u,v,pth,uh_in,sh_in) =  ((p->Em(p))∘pth) * (ε(v) ⊙ (σfun∘(ε(u), ε(uh_in), sh_in)))
a_Disp(u,v) =∫(A_Disp(u,v,pth,uh_in,sh_in))fem_params.dΩ
op_Disp = AffineFEOperator(a_Disp ,b_Disp ,fem_params.U0_Disp ,fem_params.V0_Disp)
uh_out = solve(op_Disp)
return  get_free_dof_values(uh_out)
end
vApp = 0.0
delv = 1.5e-2
const vAppMax = 1.0
vAppmid = 5e-1
delvmid = 1.5e-2
co = Int(round(((vAppmid - vApp)/ delv) + ((vAppMax - vAppmid)/ delvmid)))
innerMax = 30
count = 0
ψPlusPrev = CellState(0.0,dΩ)
p0 = ones(fem_params.np)
pf_vec = Filter(p0;r, fem_params)
pfh = FEFunction(fem_params.Pf, pf_vec)
pth = (pf -> Threshold(pf; β, η)) ∘ pfh
sPrev = CellState(1.0,dΩ)
sh = project(sPrev, model ,dΩ,order)

while vApp .< vAppMax
count = count .+ 1
if vApp >= vAppmid
delv = delvmid
end
vApp = vApp .+ delv
print("\n step", float(vApp))
for  inner = 1: innerMax
ψhPlusPrev = project(ψPlusPrev ,model ,dΩ,order)
RelErr = abs(sum(∫(Gc*ls*∇(sh)⋅∇(sh) + 2*ψhPlusPrev*sh*sh + (Gc/ls)*sh*sh)fem_params.dΩ-∫( (Gc/ls)*sh)fem_params.dΩ))/abs(sum(∫( (Gc/ls)*sh)fem_params.dΩ))
sh_vec = stepPhaseField(fem_params,uh,ψhPlusPrev)
sh = FEFunction(fem_params.U_PF, sh_vec)
u_vec = stepDisp(fem_params,pth, uh, sh, vApp)
uh = FEFunction(fem_params.U0_Disp, u_vec)
ψhPos_in = ψPos∘(ε(uh))
update_state!( new_EnergyState ,ψPlusPrev ,ψhPos_in)
if   RelErr  < 1e-8
break
end
end
end
writevtk(Ω,"resultstph",cellfields=["uh"=>uh, "sh"=>sh])
``````