Implement the dynamic J integral for phase field

[lyyc199586]

Reason

We are looking at the energy release rate in dynamic fracture now, and the dynamic J integral will be very crucial.

Design

The formula for dynamic J integral looks clear, just adding the kinetic part compared to the static one:
quasi-static:

$$ J = \int_\Gamma(Wn_1-t_i u_{i,1}) d\Gamma $$

dynamic:

$$ J = \int_\Gamma[ (W + \frac 1 2 \rho \dot{u}^2) n_1-t_i u_{i,1}] d\Gamma $$

Some problems on the implementation

Currently, I am just adding the density in the input params and add the dynamic energy density in the calculation:

#include "DynamicPhaseFieldJIntegral.h"

registerMooseObject("raccoonApp", DynamicPhaseFieldJIntegral);

InputParameters
DynamicPhaseFieldJIntegral::validParams()
{
  InputParameters params = SideIntegralPostprocessor::validParams();
  params += BaseNameInterface::validParams();
  params.addClassDescription("Compute the J integral for a phase-field model of fracture");
  params.addRequiredParam<RealVectorValue>("J_direction", "direction of J integral");
  params.addParam<MaterialPropertyName>("strain_energy_density",
                                        "psie"
                                        "Name of the strain energy density");
  params.addRequiredCoupledVar(
      "displacements",
      "The displacements appropriate for the simulation geometry and coordinate system");
  params.addParam<MaterialPropertyName>(
      "density", "density", "Name of material property containing density");
  return params;
}

DynamicPhaseFieldJIntegral::DynamicPhaseFieldJIntegral(const InputParameters & parameters)
  : SideIntegralPostprocessor(parameters),
    BaseNameInterface(parameters),
    _stress(getADMaterialPropertyByName<RankTwoTensor>(prependBaseName("stress"))),
    _psie(getADMaterialProperty<Real>(prependBaseName("strain_energy_density"))),
    _ndisp(coupledComponents("displacements")),
    _grad_disp(coupledGradients("displacements")),
    _t(getParam<RealVectorValue>("J_direction")),
    _rho(getADMaterialProperty<Real>(prependBaseName("density", true))),
    _u_dots(coupledDots("displacements"))
{
  // set unused dimensions to zero
  for (unsigned i = _ndisp; i < 3; ++i) {
    _grad_disp.push_back(&_grad_zero);
    _u_dots.push_back(&_zero);
  }
}

Real
DynamicPhaseFieldJIntegral::computeQpIntegral()
{
  // grad(u) and dot(u)
  auto H = RankTwoTensor::initializeFromRows(
      (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);
  RealVectorValue u_dot((*_u_dots[0])[_qp], (*_u_dots[1])[_qp], (*_u_dots[2])[_qp]);

  // kinetic energy density
  ADReal psik = 0.5 * raw_value(_rho[_qp]) * u_dot * u_dot;

  RankTwoTensor I2(RankTwoTensor::initIdentity);
  ADRankTwoTensor Sigma = (_psie[_qp] + psik) * I2 - H.transpose() * _stress[_qp];
  RealVectorValue n = _normals[_qp];
  
  return raw_value(_t * Sigma * n);
}

The input block:

[DJint]
    type = DynamicPhaseFieldJIntegral
    J_direction = '1 0 0'
    strain_energy_density = psie
    displacements = 'disp_x disp_y'
    density = density
    boundary = 'left bottom right top'
[]

However this gave me unreasonable results, I got negative J integral values from a simple tensile test:

Here's the PF results and energies:

However the J-integral looks wrong:

Do you have any insights on what might be causing it?