Applying Body force in peridynamics module

[Narayanan1999]

Discussed in https://github.com/idaholab/moose/discussions/26984

@lindsayad @hchen139 ^{Originally posted by **Narayanan1999** March 6, 2024} 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 []

[hchen139]

There are a couple of issues in your input file. First, the BodyForce kernel for FEA cannot directly used for PD simulation. My understanding is that FEA involves a volume integration over an element to determine the equivalent nodal force for the body force. The PD mesh only has line elements and other information such as volume of each PD node is not used in the BodeForce or its parent class to determine the equivalent nodal force. Instead, for static problem, you can use the ConstantRate nodal kernel to apply the force due to traction, where you will manually calculate the force to be applied to PD nodes. Second, for use of fictitious layer in PD, my understanding is that you only do this for Dirichlet-type BC not Neumann-type BC. Additional comment for the mesh, the current mesh generation is too complicated for this simple problem. It can be greatly simplified.

[lindsayad]

@Narayanan1999 please don't create an issue until it's clear it should be one from the Discussion. We should have had this discussion on #26984 instead of here