Why the heat released by CoupledForce and calculate by SideDiffusiveFluxIntegral are slightly different? #30060
-
Check these boxes if you have followed the posting rules.
QuestionI find it seems that the energy conservation is not precious in MOOSE? Let's consider a cylinder with cosine heat source along the axial. The size of this cylinder: radius: 0.5; height: 10.0. The heat conduction coefficient: 200. The power distribution: (2240000)(cos(pi/5(z-5))+1). Outer surface temperature: 300. So I put my input file like: I put the input file here: But the result is |
Beta Was this translation helpful? Give feedback.
Replies: 3 comments 15 replies
-
Both of these integrals are computed using a quadrature. And unlike in finite volume where the integral in the postprocessing is the exact same as the one when forming the residual, in finite element both integrals are distinct. The divergence theorem which guarantees that our volumetric integral also ensures the conservation of sum_of_surface_fluxes + volume_source/sinks suffers from this quadrature integration error. You could raise the order of the numerical quadrature in the Executioner/quadrature block. Refining further can help, though there might be a limit if there is an error from the quadrature
how big is it? |
Beta Was this translation helpful? Give feedback.
-
|
Oh! Your idea is from this page? (Executioner/Quadrature/custom_orders)(size: 2) are of different lengths I do not know something wrong with my input file.... |
Beta Was this translation helpful? Give feedback.
-
|
Rather than custom order, set the element order and side order parameters |
Beta Was this translation helpful? Give feedback.
-
|
Tell you an unfortunate story, sir: I try But the result is the same... |
Beta Was this translation helpful? Give feedback.
-
|
can you plot this error as a function of the mesh size? |
Beta Was this translation helpful? Give feedback.
-
|
You're right sir. |
Beta Was this translation helpful? Give feedback.
-
|
Oh thanks, I remember your description. But for LAGRANGE on geometric element PRISM6, the max order is FIRST, and it cannot be SECOND. For HERMITE on geometric element PRISM6, the max order is CONSTANT. And for MONOMIAL, it will edit too much time for me because the temp will become a element variable. I will try to modify it if I have time. |
Beta Was this translation helpful? Give feedback.
-
|
What you are trying to change is the order of the finite element families. This could help, but I would like you to change the order of the quadrature first |
Beta Was this translation helpful? Give feedback.
-
|
I insert in the input file |
Beta Was this translation helpful? Give feedback.
-
|
This is for second order meshes, not changing the quadrature This is for the quadrature |
Beta Was this translation helpful? Give feedback.
-
|
Dirichlet BC does not provide the exact conservation with the way how you evaluate the heat flux on the Dirichlet boundary. But it does have the consistency, i.e. when refining the mesh, the difference will be getting smaller. If conservation is a concern, you can use Robin BC. |
Beta Was this translation helpful? Give feedback.
-
|
Oh? Actually not... I change the [rightbc] into this but it does not work... eivp ∭qv dV is still larger than output_heat ∯ k∇T ndS. |
Beta Was this translation helpful? Give feedback.
-
|
The leakage on the boundary needs to be evaluated with the solution directly based on the boundary condition |
Beta Was this translation helpful? Give feedback.
-
|
You are right It gives me the conservation answer! |
Beta Was this translation helpful? Give feedback.
-
|
This is with the convective BC too right? I think we'll need to add more nuance to the SideDiffusiveFluxIntegral |
Beta Was this translation helpful? Give feedback.
-
|
Oh. We mean that there is always a slight difference between SideDiffusiveFluxIntegral (∯ k∇T ndS) and ∭ qv dV. Maybe Galerkin's method causes this? But if we use ConvectiveHeatTransferSideIntegral (∯ h(T-Tf) ndS, certainly Tf and h must be calculated from energy conservation), we can get the conservation answer. |
Beta Was this translation helpful? Give feedback.
-
|
Oh. We mean that there is always a slight difference between SideDiffusiveFluxIntegral (∯ k∇T ndS) and ∭ qv dV. Maybe Galerkin's method causes this? But if we use ConvectiveHeatTransferSideIntegral (∯ h(T-Tf) ndS, certainly Tf and h must be calculated from energy conservation), we can get the conservation answer. |
Beta Was this translation helpful? Give feedback.
The leakage on the boundary needs to be evaluated with the solution directly based on the boundary condition$h(T-T_\infty)+ k\nabla T\cdot\vec{n}=0$ , i.e. $-k\nabla T\cdot\vec{n} = h(T-T_\infty)$ . It needs to be $\oint h(T-T_\infty)dS$ but not with the gradient of the solution. If you use Dirichlet BC, an extra postprocessing step is needed to get the conservative leakage (this paper can help understanding this https://doi.org/10.1006/jcph.2000.6577).