when calculating dislocation stress at grain boundaries, densities are now extrapolated similarly to like it was already done at free surfaces

This commit is contained in:
Christoph Kords 2010-06-07 16:01:37 +00:00
parent 3989cc2688
commit 3ce1882d29
3 changed files with 25 additions and 17 deletions

View File

@ -310,7 +310,7 @@ return
endfunction endfunction
subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,ipc,ip,el) subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,disorientation,ipc,ip,el)
!********************************************************************* !*********************************************************************
!* This function calculates from state needed variables * !* This function calculates from state needed variables *
!* INPUT: * !* INPUT: *
@ -325,7 +325,8 @@ subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,ipc,ip,el)
material_phase, & material_phase, &
homogenization_maxNgrains homogenization_maxNgrains
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
mesh_maxNips mesh_maxNips, &
mesh_maxNipNeighbors
use constitutive_j2 use constitutive_j2
use constitutive_phenopowerlaw use constitutive_phenopowerlaw
use constitutive_dislotwin use constitutive_dislotwin
@ -337,6 +338,7 @@ integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: Temperature real(pReal), intent(in) :: Temperature
real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp
real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: disorientation
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
@ -350,7 +352,7 @@ real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)
call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el) call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Tstar_v, Fe, Fp, ipc, ip, el) call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Tstar_v, Fe, Fp, disorientation, ipc, ip, el)
end select end select
@ -405,7 +407,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ip
endsubroutine endsubroutine
subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fe, Fp, Temperature, misorientation, subdt, ipc, ip, el) subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fe, Fp, Temperature, disorientation, subdt, ipc, ip, el)
!********************************************************************* !*********************************************************************
!* This subroutine contains the constitutive equation for * !* This subroutine contains the constitutive equation for *
!* calculating the rate of change of microstructure * !* calculating the rate of change of microstructure *
@ -439,7 +441,7 @@ integer(pInt), intent(in) :: ipc, ip, el
real(pReal), intent(in) :: Temperature, & real(pReal), intent(in) :: Temperature, &
subdt subdt
real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: & real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: &
misorientation disorientation
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe, & Fe, &
Fp Fp
@ -466,7 +468,7 @@ select case (phase_constitution(material_phase(ipc,ip,el)))
constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, misorientation, subdt, & call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, disorientation, subdt, &
constitutive_state, constitutive_subState0, constitutive_relevantState, subdt, ipc, ip, el) constitutive_state, constitutive_subState0, constitutive_relevantState, subdt, ipc, ip, el)
end select end select

View File

@ -749,7 +749,7 @@ endfunction
!********************************************************************* !*********************************************************************
!* calculates quantities characterizing the microstructure * !* calculates quantities characterizing the microstructure *
!********************************************************************* !*********************************************************************
subroutine constitutive_nonlocal_microstructure(state, Temperature, Tstar_v, Fe, Fp, g, ip, el) subroutine constitutive_nonlocal_microstructure(state, Temperature, Tstar_v, Fe, Fp, disorientation, g, ip, el)
use prec, only: pReal, & use prec, only: pReal, &
pInt, & pInt, &
@ -768,6 +768,7 @@ use debug, only: debugger, &
selectiveDebugger selectiveDebugger
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
mesh_maxNips, & mesh_maxNips, &
mesh_maxNipNeighbors, &
mesh_element, & mesh_element, &
FE_NipNeighbors, & FE_NipNeighbors, &
mesh_ipNeighborhood, & mesh_ipNeighborhood, &
@ -797,6 +798,8 @@ real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)
Fp ! plastic deformation gradient Fp ! plastic deformation gradient
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: &
disorientation ! crystal disorientation between me and my neighbor as quaternion
!*** input/output variables !*** input/output variables
@ -831,6 +834,7 @@ real(pReal) gb, & ! short notation f
L, & ! dislocation segment length L, & ! dislocation segment length
r1_2, & r1_2, &
r2_2 r2_2
real(pReal), dimension(6) :: transmissivity ! transmissivity factor for each interface
real(pReal), dimension(2) :: lambda ! distance of (x y z) from the segment end projected onto the dislocation segment real(pReal), dimension(2) :: lambda ! distance of (x y z) from the segment end projected onto the dislocation segment
real(pReal), dimension(3,6) :: connectingVector ! vector connecting the centers of gravity of me and my neigbor (for each neighbor) real(pReal), dimension(3,6) :: connectingVector ! vector connecting the centers of gravity of me and my neigbor (for each neighbor)
real(pReal), dimension(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress in Mandel notation real(pReal), dimension(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress in Mandel notation
@ -903,9 +907,11 @@ invFe = math_inv3x3(Fe)
do n = 1,FE_NipNeighbors(mesh_element(2,el)) do n = 1,FE_NipNeighbors(mesh_element(2,el))
transmissivity(n) = constitutive_nonlocal_transmissivity(disorientation(:,n))
neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
if ( neighboring_ip == 0 ) & if ( neighboring_ip == 0 .or. transmissivity(n) < 1.0_pReal ) & ! if no neighbor present or at grain boundary, don't calculate anything, since we use mirrored connecting vector of opposite neighbor
cycle cycle
neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el)) neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el))
@ -919,8 +925,8 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
opposite_n = n - 1_pInt + 2_pInt*mod(n,2_pInt) opposite_n = n - 1_pInt + 2_pInt*mod(n,2_pInt)
opposite_el = mesh_ipNeighborhood(1,opposite_n,ip,el) opposite_el = mesh_ipNeighborhood(1,opposite_n,ip,el)
opposite_ip = mesh_ipNeighborhood(2,opposite_n,ip,el) opposite_ip = mesh_ipNeighborhood(2,opposite_n,ip,el)
if ( opposite_ip == 0 ) & ! if no opposite neighbor present... if ( opposite_ip == 0 .or. transmissivity(opposite_n) < 1.0_pReal ) & ! if no opposite neighbor present or at grain boundary ...
connectingVector(:,opposite_n) = -connectingVector(:,n) ! ... assume "ghost" IP at mirrored position of this neighbor connectingVector(:,opposite_n) = -connectingVector(:,n) ! ... use mirrored connecting vector of opposite neighbor
enddo enddo
@ -929,17 +935,17 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
if ( neighboring_ip == 0 ) then ! if no neighbor present... if ( neighboring_ip == 0 .or. transmissivity(n) < 1.0_pReal ) then ! if no neighbor present or at grain boundary ...
opposite_n = n - 1_pInt + 2_pInt*mod(n,2_pInt) opposite_n = n - 1_pInt + 2_pInt*mod(n,2_pInt)
opposite_el = mesh_ipNeighborhood(1,opposite_n,ip,el) opposite_el = mesh_ipNeighborhood(1,opposite_n,ip,el)
opposite_ip = mesh_ipNeighborhood(2,opposite_n,ip,el) opposite_ip = mesh_ipNeighborhood(2,opposite_n,ip,el)
if ( opposite_ip == 0 ) & if ( opposite_ip == 0 .or. transmissivity(opposite_n) < 1.0_pReal ) & ! (special case if no valid neighbor on both sides)
cycle cycle
neighboring_el = opposite_el neighboring_el = opposite_el
neighboring_ip = opposite_ip neighboring_ip = opposite_ip
forall (t = 1:8, s = 1:ns) & forall (t = 1:8, s = 1:ns) &
neighboring_rhoSgl(s,t) = max(0.0_pReal, & neighboring_rhoSgl(s,t) = max(0.0_pReal, &
2.0_pReal * state(g,ip,el)%p((t-1)*ns+s) - state(g,opposite_ip,opposite_el)%p((t-1)*ns+s) ) ! ... extrapolate density from opposite neighbor (but assure positive value for density) 2.0_pReal * state(g,ip,el)%p((t-1)*ns+s) - state(g,opposite_ip,opposite_el)%p((t-1)*ns+s) ) ! ... extrapolate density from opposite neighbor (but assure positive value for density)
else else
forall (t = 1:8) neighboring_rhoSgl(:,t) = state(g,neighboring_ip,neighboring_el)%p((t-1)*ns+1:t*ns) forall (t = 1:8) neighboring_rhoSgl(:,t) = state(g,neighboring_ip,neighboring_el)%p((t-1)*ns+1:t*ns)
endif endif
@ -1032,7 +1038,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
if ( selectiveDebugger .and. verboseDebugger) then if ( selectiveDebugger .and. verboseDebugger) then
write(6,*) write(6,*)
write(6,'(a20,i1,x,i2,x,i5))') '::: microstructure',g,ip,el write(6,'(a20,i1,x,i2,x,i5)') '::: microstructure ',g,ip,el
write(6,'(i2)') n write(6,'(i2)') n
write(6,'(2(a20,x,e12.3,5x))') 'delta_rho_edge:', neighboring_rhoEdgeExcess(s), 'delta_rho_screw:', neighboring_rhoScrewExcess(s) write(6,'(2(a20,x,e12.3,5x))') 'delta_rho_edge:', neighboring_rhoEdgeExcess(s), 'delta_rho_screw:', neighboring_rhoScrewExcess(s)
write(6,'(2(a20,x,f12.3,5x))') 'Nedge:', neighboring_Nedge(s), 'Nscrew:', neighboring_Nscrew(s) write(6,'(2(a20,x,f12.3,5x))') 'Nedge:', neighboring_Nedge(s), 'Nscrew:', neighboring_Nscrew(s)
@ -1331,7 +1337,7 @@ real(pReal), intent(in) :: Temperature, & ! temperat
real(pReal), dimension(6), intent(in) :: Tstar_v, & ! current 2nd Piola-Kirchhoff stress in Mandel notation real(pReal), dimension(6), intent(in) :: Tstar_v, & ! current 2nd Piola-Kirchhoff stress in Mandel notation
previousTstar_v ! previous 2nd Piola-Kirchhoff stress in Mandel notation previousTstar_v ! previous 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: & real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: &
disorientation ! crystal disorientation between me and my neighbor (axis, angle pair) disorientation ! crystal disorientation between me and my neighbor as quaternion
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe, & ! elastic deformation gradient Fe, & ! elastic deformation gradient
Fp ! plastic deformation gradient Fp ! plastic deformation gradient

View File

@ -605,7 +605,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
if (crystallite_todo(g,i,e)) then ! all undone crystallites if (crystallite_todo(g,i,e)) then ! all undone crystallites
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, &
crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states crystallite_Fp, crystallite_disorientation(:,:,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
constitutive_previousDotState2(g,i,e)%p = 0.0_pReal constitutive_previousDotState2(g,i,e)%p = 0.0_pReal
constitutive_previousDotState(g,i,e)%p = 0.0_pReal constitutive_previousDotState(g,i,e)%p = 0.0_pReal
constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotStates constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotStates
@ -1143,7 +1143,7 @@ endsubroutine
! update the microstructure ! update the microstructure
constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(:,g,i,e), crystallite_Fe, crystallite_Fp, &
g, i, e) crystallite_disorientation(:,:,g,i,e), g, i, e)
! setting flag to true if state is below relative tolerance, otherwise set it to false ! setting flag to true if state is below relative tolerance, otherwise set it to false