Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add DCDD shock capture stabilization term to heat transfer #1141

Merged
merged 48 commits into from
May 23, 2024

Conversation

voferreira
Copy link
Collaborator

@voferreira voferreira commented May 15, 2024

Description of the problem

  • Cross-wind flow leads to numerical instability in the heat transfer solver.
  • Some of this instability was not being sufficiently dampened by the GGLS term we had before, especially at highly predominant convective problems (high Pe).

Description of the solution

  • Apply the Discontinuity-Capturing Directional Dissipation (DCDD) stabilization term by Tezduyar (2003).
  • The term identifies crosswind regions and adds to them an artificial diffusivity, dampening the instabilities.

How Has This Been Tested?

  • Canonical circular advection-diffusion problem in steady state.
  • Some tests needed to be updated.

Documentation

  • To be improved

Future changes

  • Other terms exist that can be added, such as LSIC. For now, I believe this one will suffice.

Copy link
Collaborator

@hepap hepap left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice addition!

I'm just not sure I understand how the equation 79 of the paper was converted for the temperature. Maybe you can explain it to me?

source/solvers/heat_transfer_assemblers.cc Outdated Show resolved Hide resolved
Comment on lines 84 to 86
const double vdcdd = 0.5 * h * h * u_mag *
previous_temperature_gradient.norm() /
scratch_data.global_T_mag;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Which equation of the paper are you implementing here? Is it 79?


// We remove the diffusion aligned with the velocity
// as is done in the original article. In Tezduyar 2003, this is
// denoted s.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
// denoted s.
// denoted s in equation (67).

Comment on lines 93 to 95
const Tensor<2, dim> k_corr =
Utilities::fixed_power<2>(gradient_unit_vector * velocity_unit_vector) *
outer_product(velocity_unit_vector, velocity_unit_vector);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
const Tensor<2, dim> k_corr =
Utilities::fixed_power<2>(gradient_unit_vector * velocity_unit_vector) *
outer_product(velocity_unit_vector, velocity_unit_vector);
// Denoted (r.s)^2*ss in equation (70) of Tezduyar 2003.
const Tensor<2, dim> k_corr =
Utilities::fixed_power<2>(gradient_unit_vector * velocity_unit_vector) *
outer_product(velocity_unit_vector, velocity_unit_vector);

const Tensor<2, dim> k_corr =
Utilities::fixed_power<2>(gradient_unit_vector * velocity_unit_vector) *
outer_product(velocity_unit_vector, velocity_unit_vector);
const Tensor<2, dim> gradient_unit_tensor =
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
const Tensor<2, dim> gradient_unit_tensor =
// Denoted rr in equation (70) of Tezduyar 2003.
const Tensor<2, dim> gradient_unit_tensor =

outer_product(velocity_unit_vector, velocity_unit_vector);
const Tensor<2, dim> gradient_unit_tensor =
outer_product(gradient_unit_vector, gradient_unit_vector);
const Tensor<2, dim> dcdd_factor = gradient_unit_tensor - k_corr;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
const Tensor<2, dim> dcdd_factor = gradient_unit_tensor - k_corr;
// Denoted rr - (r.s)^2*ss in equation (70) of Tezduyar 2003.
const Tensor<2, dim> dcdd_factor = gradient_unit_tensor - k_corr;

// solver.
const double tolerance = 1e-12;

// In Tezduyar 2003, this is denoted r
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add the same comment as the previous for the equations and terms?

@voferreira voferreira added WIP When a PR is open but not ready for review Only documentation missing and removed Ready for review labels May 16, 2024
@voferreira
Copy link
Collaborator Author

@blaisb My bad. I will rebase from the next time on.

@blaisb
Copy link
Contributor

blaisb commented May 22, 2024

is this ready to review or? The flags says ready to review but the tests do not pass.

@voferreira
Copy link
Collaborator Author

It is ready to be reviewed. There is only one test failing. It is a new test and it is just missing updating. I will update it when I arrive at the office.

@blaisb
Copy link
Contributor

blaisb commented May 22, 2024

@voferreira What about my comment regarding convergence? I did not see any answer for the high Pe case.

Copy link
Contributor

@blaisb blaisb left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm good with the PR, but there are two points remaining
1- Why was the vof test changed to be without VOF now? If this test is without VOF now, the name of the test should be changed and the reason for this change should be explained in the PR ( I am unsure why such a change would be necessary)
2- What happens to the convergence of the high Pe case?

Pressure drop: 0 Pa
Total pressure drop: 0 Pa

----------------------
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is the test case without VOF now? I don't understand, this used to be VOF test case...

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It should still be a test for VOF. I guess the results for the same test but without VOF were put here by mistake :) it's almost the same file names! However, I think both should change with this modification.

Pressure drop: 0 Pa
Total pressure drop: 0 Pa

----------------------
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It should still be a test for VOF. I guess the results for the same test but without VOF were put here by mistake :) it's almost the same file names! However, I think both should change with this modification.

source/solvers/heat_transfer_assemblers.cc Outdated Show resolved Hide resolved
source/solvers/heat_transfer_assemblers.cc Outdated Show resolved Hide resolved
source/solvers/heat_transfer_assemblers.cc Outdated Show resolved Hide resolved
@voferreira
Copy link
Collaborator Author

@blaisb Sorry for the delay about the comment on the high Pe, but I wanted to confirm with Mikael before. The test was not passing when: 1- we were using the previous solution for steady and transient, and 2- when we were testing the term without the velocity component. Now the test is not affected anymore. As for the VOF test, it was changed by mistake and I am fixing it just now.

voferreira and others added 2 commits May 22, 2024 13:29
Co-authored-by: hepap <47506601+hepap@users.noreply.github.com>
@blaisb blaisb self-requested a review May 22, 2024 18:14
@blaisb
Copy link
Contributor

blaisb commented May 22, 2024

Good for me then, i'll wait for the other to approve before merging

Copy link
Collaborator

@PierreLaurentinCS PierreLaurentinCS left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very fine addition, good job, just switch the indices in the assembly and we're good

include/solvers/heat_transfer.h Outdated Show resolved Hide resolved
double delta_T_ref =
std::max(solution_maximum - solution_minimum, minimum_delta_T_ref);

// Set minimum value as 1 to prevent overshooting of diffusivity
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This comment should be put above the definition of delta_T_ref to explain why we use std::max, not after

source/solvers/heat_transfer_assemblers.cc Show resolved Hide resolved
source/solvers/heat_transfer_assemblers.cc Outdated Show resolved Hide resolved
source/solvers/heat_transfer_assemblers.cc Show resolved Hide resolved
source/solvers/heat_transfer_assemblers.cc Show resolved Hide resolved
source/solvers/heat_transfer_assemblers.cc Outdated Show resolved Hide resolved
source/solvers/heat_transfer_assemblers.cc Outdated Show resolved Hide resolved
source/solvers/heat_transfer_assemblers.cc Show resolved Hide resolved
source/solvers/heat_transfer_assemblers.cc Show resolved Hide resolved
Copy link
Collaborator

@PierreLaurentinCS PierreLaurentinCS left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Everything is fine, good job!

@blaisb blaisb merged commit 372edfb into master May 23, 2024
8 checks passed
@blaisb blaisb deleted the heat_shock_capture_dcdd branch May 23, 2024 22:47
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

7 participants