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.

Sign up for GitHub

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

Implement temperature-dependant solid domain contraints in VOF #1048

Merged
merged 5 commits into from
Mar 5, 2024
+712 −104
Merged

Implement temperature-dependant solid domain contraints in VOF #1048

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
450fd58
Implement temperature-dependant solid domain contraints in VOF
AmishgaAlphonius Feb 29, 2024
107c55d
Correct typos
AmishgaAlphonius Feb 29, 2024
69b8ecb
Rename subsection and function
AmishgaAlphonius Mar 4, 2024
6f0317c
Apply suggestions from code review
AmishgaAlphonius Mar 4, 2024
5e96906
Apply suggestions from code review
AmishgaAlphonius Mar 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions applications_tests/lethe-fluid/heat_transfer_phase_change_constrain_solid_domain.prm
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,10 +119,10 @@ subsection boundary conditions heat transfer
end

#---------------------------------------------------
# Constrain Solid Domain
# Constrain Stasis
#---------------------------------------------------

subsection constrain solid domain
subsection constrain stasis
set enable = true
set number of constraints = 1
subsection constraint 0
Expand Down
62 changes: 62 additions & 0 deletions applications_tests/lethe-fluid/heat_transfer_vof_phase_change_constrain_solid_domain.output
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
Running on 1 MPI rank(s)...
Number of active cells: 256
Number of degrees of freedom: 867
Volume of triangulation: 0.0625
Number of thermal degrees of freedom: 289
Number of VOF degrees of freedom: 289
Pressure drop: 0 Pa
Total pressure drop: 0 Pa

---------------
VOF Barycenter
---------------
time x_vof y_vof vx_vof vy_vof
0.000000 0.125000 0.058681 0.000000 0.000000

*******************************************************************************
Transient iteration: 1 Time: 0.1 Time step: 0.1 CFL: 0
*******************************************************************************
Pressure drop: -0.366978 Pa
Total pressure drop: -0.366978 Pa

---------------
VOF Barycenter
---------------
time x_vof y_vof vx_vof vy_vof
0.100000 0.125000 0.058681 -0.000009 -0.000023
Copy link
Collaborator

Choose a reason for hiding this comment

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

Maybe you can add also the average velocity? I think the test is still good, but maybe velocity could be a good metric too.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I tried it before, but the average velocity doesn't print anything .-.


*********************************************************************************
Transient iteration: 2 Time: 0.2 Time step: 0.1 CFL: 0.00412505
*********************************************************************************
Pressure drop: -0.367039 Pa
Total pressure drop: -0.367039 Pa

---------------
VOF Barycenter
---------------
time x_vof y_vof vx_vof vy_vof
0.200000 0.124998 0.058682 -0.000002 0.000006

*********************************************************************************
Transient iteration: 3 Time: 0.3 Time step: 0.1 CFL: 0.00230552
*********************************************************************************
Pressure drop: -0.367051 Pa
Total pressure drop: -0.367051 Pa

---------------
VOF Barycenter
---------------
time x_vof y_vof vx_vof vy_vof
0.300000 0.124995 0.058684 0.000003 0.000031

********************************************************************************
Transient iteration: 4 Time: 0.4 Time step: 0.1 CFL: 0.0016325
********************************************************************************
Pressure drop: -0.367145 Pa
Total pressure drop: -0.367145 Pa

---------------
VOF Barycenter
---------------
time x_vof y_vof vx_vof vy_vof
0.400000 0.124992 0.058687 0.000002 0.000032
231 changes: 231 additions & 0 deletions applications_tests/lethe-fluid/heat_transfer_vof_phase_change_constrain_solid_domain.prm
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,231 @@
# Listing of Parameters
#----------------------

set dimension = 2

#---------------------------------------------------
# Simulation Control
#---------------------------------------------------

subsection simulation control
set method = bdf1
set time end = 0.4
set time step = 0.1
set output frequency = 1
AmishgaAlphonius marked this conversation as resolved.
Show resolved Hide resolved
end

#---------------------------------------------------
# Multiphysics
#---------------------------------------------------

subsection multiphysics
set heat transfer = true
set fluid dynamics = true
set VOF = true
end

#---------------------------------------------------
# Initial condition
#---------------------------------------------------

subsection initial conditions
subsection temperature
set Function expression = (40-25)*(0.25-x)*4+25
end
subsection VOF
set Function expression = if(y>=0.125, 0, 1)
end
end

#---------------------------------------------------
# Source term
#---------------------------------------------------

subsection source term
subsection xyz
set Function expression = 5 ; 0 ; 0
end
end

#---------------------------------------------------
# Physical Properties
#---------------------------------------------------

subsection physical properties
set number of fluids = 2
set reference temperature = 30
subsection fluid 0
set thermal conductivity = 0.5
set thermal expansion model = phase_change
set rheological model = phase_change
set specific heat model = phase_change
subsection phase change
set latent enthalpy = 200
set liquidus temperature = 30
set solidus temperature = 29.8
set viscosity liquid = 0.0001
set viscosity solid = 10
end
end
subsection fluid 1
set thermal conductivity = 0.5
set thermal expansion model = phase_change
set rheological model = phase_change
set specific heat model = phase_change
subsection phase change
set latent enthalpy = 200
set liquidus temperature = 35
set solidus temperature = 34.8
set viscosity liquid = 0.0001
set viscosity solid = 10
end
end
end

#---------------------------------------------------
# Mesh
#---------------------------------------------------

subsection mesh
set type = dealii
set grid type = hyper_cube
set grid arguments = 0 : 0.25 : true
set initial refinement = 4
end

#---------------------------------------------------
# Boundary Conditions
#---------------------------------------------------

subsection boundary conditions
set number = 4
subsection bc 0
set id = 0
set type = noslip
end
subsection bc 1
set id = 1
set type = noslip
end
subsection bc 2
set id = 2
set type = noslip
end
subsection bc 3
set id = 3
set type = noslip
end
end

subsection boundary conditions heat transfer
set number = 2
subsection bc 0
set id = 0
set type = temperature
subsection value
set Function expression = 40
end
end
subsection bc 1
set id = 1
set type = temperature
subsection value
set Function expression = 25
end
end
end

#---------------------------------------------------
# Constrain Stasis
#---------------------------------------------------

subsection constrain stasis
set enable = true
set number of constraints = 2
subsection constraint 0
set fluid id = 0
set phase fraction tolerance = 1e-4
set min temperature = 0
set max temperature = 29
end
subsection constraint 1
set fluid id = 1
set phase fraction tolerance = 1e-4
set min temperature = 0
set max temperature = 34
end
end

#---------------------------------------------------
# Post-processing
#---------------------------------------------------

subsection post-processing
set verbosity = verbose
set calculate pressure drop = true
set calculate mass conservation = false
AmishgaAlphonius marked this conversation as resolved.
Show resolved Hide resolved
set calculate barycenter = true
end

#---------------------------------------------------
# Non-Linear Solver Control
#---------------------------------------------------

subsection non-linear solver
subsection fluid dynamics
set tolerance = 1e-4
set max iterations = 5
set verbosity = quiet
end
subsection heat transfer
set tolerance = 1e-3
set max iterations = 20
set verbosity = quiet
end
subsection VOF
set tolerance = 1e-9
set max iterations = 20
set verbosity = quiet
end
end

#---------------------------------------------------
# Linear Solver Control
#---------------------------------------------------

subsection linear solver
subsection fluid dynamics
set verbosity = quiet
set method = gmres
set max iters = 1000
set relative residual = 1e-3
set minimum residual = 1e-5
set preconditioner = ilu
set ilu preconditioner fill = 0
set ilu preconditioner absolute tolerance = 1e-12
set ilu preconditioner relative tolerance = 1.00
set max krylov vectors = 1000
end
subsection heat transfer
set verbosity = quiet
set method = gmres
set max iters = 1000
set relative residual = 1e-3
set minimum residual = 1e-5
set preconditioner = ilu
set ilu preconditioner fill = 0
set ilu preconditioner absolute tolerance = 1e-12
set ilu preconditioner relative tolerance = 1.00
end
subsection VOF
set verbosity = quiet
set method = gmres
set max iters = 1000
set relative residual = 1e-3
set minimum residual = 1e-10
set preconditioner = ilu
set ilu preconditioner fill = 0
set ilu preconditioner absolute tolerance = 1e-12
set ilu preconditioner relative tolerance = 1.00
end
end
2 changes: 1 addition & 1 deletion ...w/3d-mixer-using-single-rotating-frame/3d-mixer-using-single-rotating-frame.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ The rotating Lagrangian frame of reference is non-Galilean. Consequently, the Co
\nabla \cdot \mathbf{u} &= 0 \\
\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} &= -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u} +\mathbf{f} - \underbrace{\Omega \times \mathbf{u}}_{Coriolis} - \underbrace{\Omega \times (\mathbf{q} \times \mathbf{u})}_{Centrifugal}

where :math:`\mathbf{q}` is the position in the fluid with respect to the center of rotation and :math:`\mathbf{\Omega}` is the angular velocity of the rotating reference frame. The Coriolis force adds a velocity dependant force to the Navier-Stokes equations whereas the centrifugal forces is independent of the flow and only modifies the pressure field.
where :math:`\mathbf{q}` is the position in the fluid with respect to the center of rotation and :math:`\mathbf{\Omega}` is the angular velocity of the rotating reference frame. The Coriolis force adds a velocity dependent force to the Navier-Stokes equations whereas the centrifugal forces is independent of the flow and only modifies the pressure field.

In this example, we will start by simulating the case when :math:`Re = 1` and then follow with simulations for :math:`Re` values ranging from :math:`0.1` to :math:`100` to generate :math:`N_p` vs :math:`Re` curves.

Expand Down
2 changes: 1 addition & 1 deletion ...on-in-a-cylindrical-vessel/photon-count-calculation-in-a-cylindrical-vessel.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -243,7 +243,7 @@ The Case II figure shows the evolution of the photon count in absence of attenua

The Case III figure depicts the evolution of the photon count in absence of the attenuation due to the medium found inside the vessel and the vessel's wall. Since we use the same set of positions in all cases, :math:`\omega(\alpha)` and :math:`\omega(\theta)` remain the same for each given position of the tracer particle. The length of the path traveled by the photon inside the detector should also be the same since the same seed number is used. As seen on the Case III figure, when the particle is aligned with the axis of symmetry of the detector, the photon count reaches a maximum. At that position, the evolution of the product :math:`\omega(\alpha) \cdot \omega(\theta)` seen on the Case II figure also reaches a maximum. And the distance :math:`d(\alpha,\theta)` reaches a local maximum at that position. On the case III figure, we notice that the inflection points at :math:`y \approx -5.5` cm and at :math:`y \approx -3.7` cm (not too far from the edge of the detector's face), seen on the Case II figure, are not present anymore. This means that when :math:`y \in ]-10, -3.8[` cm, when the particle sees both the face and the lateral sides of the detector and as the particle approaches the detector's face, the distance :math:`d(\alpha,\theta)` increases making the count increase. And when :math:`y \in ]-3.8, -1.5[` cm the distance :math:`d(\alpha,\theta)` decreases in such way that it counters the rapid increase in weighting factors giving the evolution of the photon count a more parabolic shape. Finally, between :math:`y \in ]-1.5, 0]` cm, :math:`d(\alpha,\theta)` increases until reaching a local maximum.

The last case studied (Case IV) shows the evolution of the photon count when :math:`\mu_d` is so great that :math:`f_d(\alpha_j, \theta_j)` tends to :math:`1 \ \forall y \in ]-10, 10[` cm. By doing so, we can see the evolution of the count when the efficiency is independent of the interaction between the emitted :math:`\gamma`-ray and the detector. With this case, we isolate the effect of the evolution of :math:`f_a(\alpha, \theta)` on the count. More specifically, we're looking at the evolution of :math:`e(\alpha,\theta)` as the particle travels in the vessel, since :math:`\mu_r` remains constant in the studied domain. We notice that we have a local minimum at :math:`y \approx -4.6` where we saw the convex section on the Case II figure. Considering the Case II results, we can interpret the Case IV figure as follows. Starting from the back of the vessel, where :math:`f_a(\alpha, \theta)` is at its maximal value, :math:`f_a(\alpha, \theta)` decreases at a decreasing rate until reaching :math:`y \approx -4.6` cm. The maximal value of :math:`f_a(\alpha, \theta)` (minimal value of :math:`e(\alpha,\theta)`) being when the particle is the furthest away from the detector may be explained by the curvature of the vessel's wall. Since the wall of the vessel is curved to form a circle, the distance traveled by the photon inside the vessel on the average probable path isn't necessarily larger than the radius of the reactor. We know that at :math:`y = 0`, :math:`e(\alpha,\theta) = 10` cm. In other words, :math:`e(\alpha,\theta)` is equivalent to the radius of the reactor. On the :math:`e(\alpha,\theta)` *function of* :math:`y` figure, we can read :math:`e(\alpha,\theta) \approx 10.04` cm when :math:`y = 10` cm. We also know that an increasing distance :math:`e(\alpha,\theta)` leads to a decreasing efficiency, which means a decreasing count. Therefore, we may assume that :math:`e(\alpha,\theta)` is minimal when :math:`y \approx -10` cm or when :math:`y \approx 10` cm. And, it slowly increases until reaching :math:`y \approx -4.6` cm. When the particle reaches the :math:`y \approx -4.6` cm position (local minimum), the variation of :math:`f_a(\alpha, \theta)` is so little that :math:`f_a(\alpha, \theta)` behaves as a constant. This explains why we see the same pattern of evolution of the photon count as in Case II when :math:`y \in ]-4.6, -3.8[` cm. Similarly, when the particle sees only the face of the detector, the pattern of the counts evolution follows the same trend as the one seen on the Case II when :math:`y \in ]-3.8, 0]` cm. This also indicates very little fluctuations of :math:`e(\alpha,\theta)` as we may see on the :math:`e(\alpha,\theta)` *function of* :math:`y` figure. Therefore, the photon count is highly dependant of the weighting factors when :math:`y \in ]-3.8, 0]` cm.
The last case studied (Case IV) shows the evolution of the photon count when :math:`\mu_d` is so great that :math:`f_d(\alpha_j, \theta_j)` tends to :math:`1 \ \forall y \in ]-10, 10[` cm. By doing so, we can see the evolution of the count when the efficiency is independent of the interaction between the emitted :math:`\gamma`-ray and the detector. With this case, we isolate the effect of the evolution of :math:`f_a(\alpha, \theta)` on the count. More specifically, we're looking at the evolution of :math:`e(\alpha,\theta)` as the particle travels in the vessel, since :math:`\mu_r` remains constant in the studied domain. We notice that we have a local minimum at :math:`y \approx -4.6` where we saw the convex section on the Case II figure. Considering the Case II results, we can interpret the Case IV figure as follows. Starting from the back of the vessel, where :math:`f_a(\alpha, \theta)` is at its maximal value, :math:`f_a(\alpha, \theta)` decreases at a decreasing rate until reaching :math:`y \approx -4.6` cm. The maximal value of :math:`f_a(\alpha, \theta)` (minimal value of :math:`e(\alpha,\theta)`) being when the particle is the furthest away from the detector may be explained by the curvature of the vessel's wall. Since the wall of the vessel is curved to form a circle, the distance traveled by the photon inside the vessel on the average probable path isn't necessarily larger than the radius of the reactor. We know that at :math:`y = 0`, :math:`e(\alpha,\theta) = 10` cm. In other words, :math:`e(\alpha,\theta)` is equivalent to the radius of the reactor. On the :math:`e(\alpha,\theta)` *function of* :math:`y` figure, we can read :math:`e(\alpha,\theta) \approx 10.04` cm when :math:`y = 10` cm. We also know that an increasing distance :math:`e(\alpha,\theta)` leads to a decreasing efficiency, which means a decreasing count. Therefore, we may assume that :math:`e(\alpha,\theta)` is minimal when :math:`y \approx -10` cm or when :math:`y \approx 10` cm. And, it slowly increases until reaching :math:`y \approx -4.6` cm. When the particle reaches the :math:`y \approx -4.6` cm position (local minimum), the variation of :math:`f_a(\alpha, \theta)` is so little that :math:`f_a(\alpha, \theta)` behaves as a constant. This explains why we see the same pattern of evolution of the photon count as in Case II when :math:`y \in ]-4.6, -3.8[` cm. Similarly, when the particle sees only the face of the detector, the pattern of the counts evolution follows the same trend as the one seen on the Case II when :math:`y \in ]-3.8, 0]` cm. This also indicates very little fluctuations of :math:`e(\alpha,\theta)` as we may see on the :math:`e(\alpha,\theta)` *function of* :math:`y` figure. Therefore, the photon count is highly dependent of the weighting factors when :math:`y \in ]-3.8, 0]` cm.

Coming back to the Case I figure, we can see that photon count follows a pattern similar to the one seen in Case IV. We may interpret from it that :math:`f_d(\alpha, \theta)` varies very little as opposed to :math:`f_a(\alpha, \theta)` that fluctuates greatly. The local minimal values, in this case, are at :math:`y \approx -6` cm and :math:`y \approx 6` cm, as opposed to :math:`y \approx -4.6` cm and :math:`y \approx -4.6` cm for the fourth case. This is due to the change in the value of :math:`\mu_d`. :math:`f_d(\alpha,\theta)` function of :math:`y` increases at a slower rate, making the minimums further way from the center. To summarize, the fluctuations seen in the Case I figure is the result of the combined influence of the values of the attenuation coefficients, the variation of the path lengths of the photon in the vessel and the detector, and the evolution of the weighting factors.

Expand Down
2 changes: 1 addition & 1 deletion doc/source/parameters/cfd/analytical_solution.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ If there is an analytical solution for the Cahn-Hilliard physics, enter the ``ca
* The ``Function expression`` parameter sets the expression of the phase order parameter and the chemical potential.

.. note::
The variables *x*, *y*, *z* (3D) and *t* (time-dependant) can be used in the function expressions.
The variables *x*, *y*, *z* (3D) and *t* (time-dependent) can be used in the function expressions.

In all four last subsections, you can add a ``Function constants`` parameter that will act as a constant in the ``Function expression``.

Expand Down
2 changes: 1 addition & 1 deletion doc/source/parameters/cfd/cfd.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ General, CFD and Multiphysics
boundary_conditions_multiphysics
box_refinement
cahn_hilliard
constrain_solid_domain
constrain_stasis
dimensionality
dynamic_flow_control
fem
Expand Down
Loading
Loading