I think, the remedy may be first to check your equations and than to check you codes step by step very carefully. Also have a look at some online examples, e.g.,

https://classes.engineering.wustl.edu/2009/spring/mase5513/abaqus/docs/v6.6/books/exa/default.htm?startat=ch01s04aex46.html

https://abaqus-docs.mit.edu/2017/English/SIMACAEBMKRefMap/simabmk-c-jintegral2d.htm

Does the value of J-integral obtained from history output type J-integral differ significantly from its value from history output type Stress intensity factor ? Normally they are very close.

Muhammad Ali, Thank you very much for your input. I appreciate your time. I learned that actually there is a way to calculate mode I and mode II stress intensity factors from the calculated stress values around the crack tip elements for 3D cases with bending load.

Victor J. Аdlucky, Thank you very much for your input. I appreciate your time.

The J integral from history output is reasonable. But the stress intensity factors are much different from the J integral one.

There is no out of plane share stress in my case. So, I discarded K III.

I calculated the equivalent stress intensity factors from the below way in the attached PPT file.

Dear Subrato Roy,

Thanks for sharing your results. If the units of J-integral in your presentation are kN/m, then values of K correspond to the values of J-integral (I assumed that E=200GPa). Also, obtained values of K are slightly lower than critical values for steel.