Applying Body force in peridynamics module

[Narayanan1999]

Sir, I am following this book 'Peridynamic Theory and Its Applications' which suggest that the tractions are applied in terms of body force. I wanted to apply traction boundary condition on 1D bar. For that I created a bar of 10mm length. I also created a fictitious boundary layer on left end of the bar which is fixed. At the right edge, I created a separate block to apply the body force equivalent to traction boundary. I am using BPD with horizon number =1. Distance between two material point is 1mm. The resulting displacment is varying linearly as expected but it is not matching with the analytical solution. Also, if I increase the mesh, then displacement value at the end is changing in the order of around 10.
Is there a way to use Body Force kernel properly for peridynamics module? I am attaching my code along with this.

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
[]

L=10
H=1
W=1
Hor_num=1
Hn=1
dimension=3
xn=10
yn=1
zn=1
force_dV=100e4

delx=${fparse L/xn}
dely=${fparse H/yn}
delz=${fparse W/zn}
#dV=${fparse delx*dely*delz}
#force_dV= ${fparse force}
M_L_min=${fparse delx/2}
M_L_max=${fparse L-delx/2}
M_H_min=${fparse -dely/2}
M_H_max=${fparse H+dely/2}
M_W_min=${fparse -delz/2}
M_W_max=${fparse W+delz/2}
M_nx=${fparse xn-1}
fb_nx=${fparse Hn}
M_ny=${fparse yn+1}
M_nz=${fparse zn+1}



L_L_min=${fparse M_L_min-Hn*delx}
R_L_max=${fparse M_L_max+delx}

[Mesh]
  type = PeridynamicsMesh
  horizon_number = ${Hor_num}
  [Left_b]
    type = GeneratedMeshGenerator
    dim = ${dimension}
    xmin = ${L_L_min}
    xmax = ${M_L_min}
    ymin = ${M_H_min}
    ymax = ${M_H_max}
	zmin = ${M_W_min}
    zmax = ${M_W_max}
    nx = ${fb_nx}
    ny = ${M_ny}
	nz = ${M_nz}
	elem_type=HEX8
    boundary_name_prefix = Left_b
	boundary_id_offset = 10
  []
  [Left_b_id]
    type = SubdomainIDGenerator
    input = Left_b
    subdomain_id = 1
  []
  
  [Mainsec]
    type = GeneratedMeshGenerator
    dim = ${dimension}
    xmin = ${M_L_min}
    xmax = ${M_L_max}
    ymin = ${M_H_min}
    ymax = ${M_H_max}
	zmin = ${M_W_min}
    zmax = ${M_W_max}
    nx = ${M_nx}
    ny = ${M_ny}
	nz = ${M_nz}
	elem_type=HEX8
    boundary_name_prefix = Mainsec
    boundary_id_offset = 20
  []
  [Mainsec_id]
    type = SubdomainIDGenerator
    input = Mainsec
    subdomain_id = 2
  []
  [Right_b]
    type = GeneratedMeshGenerator
    dim = ${dimension}
    xmin = ${M_L_max}
    xmax = ${R_L_max}
    ymin = ${M_H_min}
    ymax = ${M_H_max}
	zmin = ${M_W_min}
    zmax = ${M_W_max}
    nx = 1
    ny = ${M_ny}
	nz = ${M_nz}
	elem_type=HEX8
    boundary_name_prefix = Right_b
	boundary_id_offset = 30
  []
  [Right_b_id]
    type = SubdomainIDGenerator
    input = Right_b
    subdomain_id = 3
  []
  [combined]
    type = MeshCollectionGenerator
    inputs = 'Left_b_id Mainsec_id Right_b_id'
  []
  construct_side_list_from_node_list=true
  [./gpd]
    type = MeshGeneratorPD
    input = combined
    retain_fe_mesh = false
	blocks_to_pd = '1 2 3'
    merge_pd_blocks = false
    bonding_block_pairs = '1 2; 2 3'
    merge_pd_interfacial_blocks = true
  [../]
[]

[Variables]
  [./disp_x]
  [../]
  [./disp_y]
  [../]
  [./disp_z]
  [../]
[]

[AuxVariables]
  [./node_vol]
  [../]
  [./s11]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[AuxKernels]
  [./node_area]
    type = NodalVolumePD
    variable = node_vol
  [../]
  [./s11]
    type = NodalRankTwoPD
    variable = s11

	youngs_modulus = 100e9
    poissons_ratio = 0

    rank_two_tensor = stress
	output_type = component
    index_i = 0
    index_j = 0
  [../]
[]

[BCs]
  [./left_x]
    type = DirichletBC
    variable = disp_x
    boundary = 998
    value = 0.0
  [../]
  [./left_y]
    type = DirichletBC
    variable = disp_y
    boundary = 998
    value = 0.0
  [../]
  [./left_z]
    type = DirichletBC
    variable = disp_z
    boundary = 998
    value = 0.0
  [../]
[]



[Kernels]
  [./force_x2]
    type = BodyForce
    variable = disp_x
    block = '1003'
    value = ${force_dV}
  [../]
[]


[Modules/Peridynamics/Mechanics/Master]
  [./all]
    formulation = BOND
  [../]
[]

[Materials]
  [./elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 100e9
    poissons_ratio = 0
  [../]

  [./force_density]
    type = ComputeSmallStrainVariableHorizonMaterialBPD
  [../]
  
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]
[Executioner]
  type = Steady
  solve_type= NEWTON
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu NONZERO'
[]

[Outputs]
  exodus = true
[]

[GiudGiud]

@hchen139

[lindsayad]

For better or worse this conversation has been forked onto #26985 where @hchen139 has replied. I'm going to close this discussion