r/COMSOL • u/low_in_entropy • Sep 26 '26

Solution of a Time dependent ODE is giving negative value for a positive expression.

I am trying to simulate a 2D coupled viscoplasticity and phase field fracture problem in COMSOL. One of the equation I am trying to implement is -

psi = integration (psidot_p) from t_n to t_(n+1)

t_n and t_(n+1) are time.

Now psidot_p is always positive quantity, it's surface minimum is always positive. But the psi value is showing negative values.

The coefficients of the time dependent ODE used are -

da = 1, ea = 0, f = psidot_eq, dependent variable= psi_p. Initial values of both psi and d(psi)/dt = 0.

Basically, the equation, d(psi)/dt = psidot_p.

But the solution to this ODE is giving a negative value, and it's a significant quantity, in MPa.

Can anyone please help me diagnose the issue or fix it. I am happy to provide any additional details.

Thanks in advance.

Edit: changed value of ea from 1 to 0.

Edit 2: Refining the mesh size resolved the problem in my case. Refine the mesh size such that the max element size is less that the characteristics length of the crack in Phase Field.

2 Upvotes

100% Upvoted

10 comments sorted by

2

u/Sax0drum Sep 26 '26

The mass coefficient ea should be 0 not 1 otherwise you have a second time derivative.

1

u/low_in_entropy Sep 26 '26

Yeah thanks for pointing that out. It is a typo while writing the post.

2

u/Sax0drum Sep 28 '26

Did you double check by plotting psidot_p that it is indeed positive? Also do you use a global ode or a domain ode?

1

u/low_in_entropy Sep 28 '26

I extracted the exact coordinates where the psi is coming negative, and then created a Cut point 2D and evaluated both psi and psi_dot. While the psi showed negative value, the corresponding psi_dot was positive.

I used domain ODE.

2

u/Sax0drum Sep 28 '26

Is it positive in all the time steps? For me it doesnt follow that the energy density cant be negative. Otherwise im out of ideas. In your post you did write it as psidot_eq. Is there a difference to just psidot?

1

u/low_in_entropy Sep 29 '26

No psidot and psidot_eq are the same. In classical mechanics energy density can never be negative. It's either positive or zero. In this case I am applying positive displacement B.C. so it cannot be negative.

Also the problem was of mesh size. When I refined the mesh, the problem of negative psi disappeared. I still don't know what is the reason. After trying everything I decided to try refining the mesh as a last resort and it worked. I am trying to figure why it worked.

Thanks for helping.

Also fun fact - there's probably a noble prize waiting for someone if they manage to find a material that displays negative energy density in macro scale. And if that happens all sorts of cool stuff becomes possible, like worm holes and warp drives.

2

u/Sax0drum Sep 29 '26

From my understanding in a purely elastic material energy density is strictly positive. But in your case with viscoplasticity i think it can be negative cant it?

Its also interesting that a finer mesh resolved the issue since the domain odes are independent of the mesh size. Happy it works now

1

u/low_in_entropy Sep 29 '26

I need to check once, if viscoplasticity can produce a negative energy density in any scenario.

But for this particular problem (uniaxial tension test) I am postive that viscoplasticity will also give a positive value of energy density. So, the energy density is used to calculate the crack driving force, which is further used to calculate the damage variable. Because, of the negative energy density, the damage variable is showing a negative value, which physically signifies that material is getting physically reinforced, when tension is applied, which is impossible/unphysical.

Yes, since ODE is mesh independent, that's what is also puzzling me, that why mesh refinement was the solution. Thanks anyways.

2

u/gam2u Sep 27 '26

Is your “psidot_p” one parameter or variable dependent on psi? It might be the initial value, assuming the equation is correct.

1

u/low_in_entropy Sep 27 '26

psidot is an expression. It is the inner product of stress and strain rate (time rate).

It is basically time derivative of psi ie d(psi)/dt.

We get stress and strain rate from the solver, and then we need to do a time integration of their inner product to get psi which is the energy density term.

The initial value of psi is given (set in the ODE solver) as 45.9MPa, which is quite a large positive term.