Issue with coupledDiv for Hydrogen Embrittlement Governing Equation

[AlejandroRiano2023]

Check these boxes if you have followed the posting rules.

• Q&A General is the most appropriate section for my question
• I have consulted the posting Guidelines on the Discussions front page
• I have searched the Discussions forum and my question has not been asked before
• I have searched the MOOSE website and the documentation does not answer my question
• I have formatted my post following the posting guidelines (screenshots as a last resort, triple back quotes around pasted text)

Question

Hello, I hope you are doing great!

Hello,

I hope you are doing well!

I am trying to implement the following governing equation for Hydrogen Embrittlement:

The equation involves the use of divergence. I am attempting to use coupledDiv as suggested in the MOOSE documentation, but I encounter the following error when compiling my application:

/home/lriaomol/projects/otter/src/kernel/KernelEq1.C:51:19: error: 'coupledDiv' was not declared in this scope; did you mean 'coupledDot'?
51 | Real term2 = -coupledDiv((1.0 - _d[_qp]) * (1.0 - _d[_qp]) * _D_L * _grad_C_L[_qp]);
| ^~~~~~~~~~
| coupledDot

Could you please assist in resolving this issue? Additionally, I am using the article “Phase Field Modeling of Hydrogen Embrittlement” by Huang (2020) as the basis for my model, which you can find here.

Here is the kernelEq1.c

#include "KernelEq1.h"
#include "MooseVariable.h"

registerMooseObject("otterApp", KernelEq1);

InputParameters
KernelEq1::validParams()
{
InputParameters params = Kernel::validParams();
params.addRequiredCoupledVar("C_L", "Concentration C_L.");
params.addRequiredCoupledVar("C_T", "Concentration C_T.");
params.addRequiredCoupledVar("theta_T", "Theta_T variable.");
params.addRequiredCoupledVar("d", "d variable.");
params.addRequiredCoupledVar("hidrostatic_stress", "Hydrostatic stress variable.");
params.addRequiredCoupledVar("N_T", "N_T variable.");
params.addRequiredCoupledVar("effective_plastic_strain", "Effective plastic strain.");
params.addParam("D_L", 12700.0, "Diffusion coefficient D_L (micrometers squared per second).");
params.addParam("V_H", 2e12, "Volume or velocity V_H (micrometers cubed per mole).");
params.addParam("R", 4.1240, "R value.");
params.addParam("Theta", 554.25, "Theta value.");
return params;
}

KernelEq1::KernelEq1(const InputParameters & parameters) :
Kernel(parameters),
_C_L(coupledValue("C_L")),
_grad_C_L(coupledGradient("C_L")),
_der_C_L(coupledDot("C_L")),
_C_T(coupledValue("C_T")),
_theta_T(coupledValue("theta_T")),
_d(coupledValue("d")),
_D_L(getParam("D_L")), // Get D_L from the parameters
_V_H(getParam("V_H")), // Get V_H from the parameters
_hidrostatic_stress(coupledValue("hidrostatic_stress")),
_grad_hidrostatic_stress(coupledGradient("hidrostatic_stress")),
_N_T(coupledValue("N_T")),
_effective_plastic_strain(coupledValue("effective_plastic_strain")),
_der_effective_plastic_strain(coupledDot("effective_plastic_strain")),
_R(getParam("R")),
_Theta(getParam("Theta"))
{
}

Real
KernelEq1::computeQpResidual()
{
// Term 1: Temporal derivative of C_L multiplied by a factor
Real term1 = ((_C_L[_qp] + _C_T[_qp] * (1.0 - _theta_T[_qp])) / _C_L[_qp]) * _der_C_L[_qp];

// Term 2: Divergence of the diffusion term
Real term2 = -coupledDiv((1.0 - _d[_qp]) * (1.0 - _d[_qp]) * _D_L * _grad_C_L[_qp]);

// Term 3: Divergence of the additional term with concentration, velocity, and gradient of sigma_H
Real term3 = coupledDiv((1.0 - _d[_qp]) * (1.0 - _d[_qp]) * _D_L * _C_L[_qp] * _V_H * _grad_hidrostatic_stress[_qp] / (_R * _Theta));

// Term 4: Temporal derivative of epsilon multiplied by the derivative of N_T with respect to epsilon
Real derivative_N_T = 2.30258 * _N_T[_qp] * 12.815 * std::exp(-5.5 * _effective_plastic_strain[_qp]);  // Derivative of N_T with respect to epsilon
Real term4 = _theta_T[_qp] * derivative_N_T * _der_effective_plastic_strain[_qp];

return term1 + term2 + term3 + term4;

}

Real
KernelEq1::computeQpJacobian()
{
// Here you can calculate the derivative with respect to the variables if necessary
return 1.0; // Simplified Jacobian
}

[GiudGiud]

Hello

coupledDiv should be used the same way that coupledValue is used, in the constructor intializer's list, and you have to specify the name of a variable, you cannot just specify an entire expression to take the divergence of

[GiudGiud]

yes exactly.
If term3 should be a vector, then use this:

RealVectorValue term3 = ....

else you are likely missing a dot-product or a "take the component of this vector" operation in the code

[AlejandroRiano2023]

Thank you so much for your response, GiudGiud! It worked perfectly.

I also have a question regarding the correct way to register the divergence of a variable when working with kernel.h. I was initially tempted to use VariableDiv, but it didn't work as expected.

Here is the relevant code:

{
public:
static InputParameters validParams();
CLCalc(const InputParameters & parameters);

protected:
/// ADKernel objects must override precomputeQpResidual
virtual ADRealVectorValue precomputeQpResidual() override;

private:
const VariableValue & _C_L; // Concentration C_L
const VariableGradient & _grad_C_L; // Gradient of C_L
const VariableDiv & _div_CLCalc_3; // Divergence of C_L (Issue here)
};

Could you help clarify the correct approach for registering the divergence of a variable?

[GiudGiud]

if you look at Coupleable.h you should see the expected type is VectorVariableDivergence

  /**
   * Returns divergence of a coupled variable
   * @param var_name Name of coupled variable
   * @param comp Component number for vector of coupled variables
   * @return Reference to a VectorVariableDivergence containing the divergence of the coupled
   * variable
   * @see Kernel::_div_u
   */
  virtual const VectorVariableDivergence & coupledDiv(const std::string & var_name,
                                                      unsigned int comp = 0) const;

[AlejandroRiano2023]

Hello GiudGiud, I reviewed the MOOSE documentation for Coupleable.h (https://mooseframework.inl.gov/releases/moose/v1.0.0/doxygen/moose/Coupleable_8h_source.html), but I couldn’t locate the expected type for VectorVariableDivergence. I don’t know if maybe I am overlooking something. Could you please share the link to the Coupleable.h file where this type is defined?

[GiudGiud]

you used the old doco from a release. Use the current one
https://mooseframework.inl.gov/docs/doxygen/moose/Coupleable_8h_source.html

[AlejandroRiano2023]

Thank you so much for your response, @GiudGiud, I really appreciate it! I just have a quick clarification: In my system of equations, I have 4 unknowns and 4 equations, and I am creating one kernel for each unknown. When I compile the app, will it automatically solve for all 4 unknowns? Is this the correct approach?

[GiudGiud]

once you form all 4 equations in the input file, whether using a single kernel per equation or multiple kernels (one per term in the equation). you will be able to run the input file using the app executable.

Note that time derivative terms need to be separated from the rest of the equation, in their own kernel

[AlejandroRiano2023]

Thank you for the clarification!

Regarding this, would it be okay to have some terms as auxiliary kernels and others as main kernels, or would that pose any issues? I'm wondering if this would affect the solution process or cause any unintended problems.

[GiudGiud]

it's ok to compute some terms with auxiliary kernels, but very often you end up needing the derivatives of these terms (either using AD or manual Jacobian). And at that point, it's better to use material properties for them

[AlejandroRiano2023]

Thank you so much for all your help, I appreciate it.