From 905fa596c8eeab6fa5409f6b4b5f44c4c721de47 Mon Sep 17 00:00:00 2001 From: Tulio Date: Mon, 11 Mar 2024 13:12:31 -0500 Subject: [PATCH] Fix thermal interface for orthotropic cases (#1009) --- .../thermally_coupled_fluid_wall.py | 28 +++-- test/test_diffusion.py | 11 +- test/test_multiphysics.py | 106 ++++++++++++++++-- 3 files changed, 116 insertions(+), 29 deletions(-) diff --git a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py index 0636c152e..392ce95ad 100644 --- a/mirgecom/multiphysics/thermally_coupled_fluid_wall.py +++ b/mirgecom/multiphysics/thermally_coupled_fluid_wall.py @@ -4,7 +4,7 @@ Couples a fluid subdomain governed by the compressible Navier-Stokes equations (:mod:`mirgecom.navierstokes`) with a wall subdomain governed by the heat equation (:mod:`mirgecom.diffusion`) through temperature and heat flux. This -radiation can optionally include a sink term representing emitted radiation. +coupling can optionally include a sink term representing emitted radiation. In the non-radiating case, coupling enforces continuity of temperature and heat flux .. math:: @@ -234,9 +234,15 @@ def viscous_divergence_flux( lengthscales_minus = project_from_base( dcoll, dd_bdry, self._lengthscales_minus) - tau = ( - self._penalty_amount * state_bc.thermal_conductivity - / lengthscales_minus) + if isinstance(state_bc.thermal_conductivity, np.ndarray): + # orthotropic materials + actx = state_minus.array_context + normal = actx.thaw(dcoll.normal(dd_bdry)) + kappa_bc = np.dot(normal, state_bc.thermal_conductivity*normal) + else: + kappa_bc = state_bc.thermal_conductivity + + tau = self._penalty_amount * kappa_bc / lengthscales_minus t_minus = state_minus.temperature t_plus = self.temperature_plus( @@ -256,7 +262,7 @@ def __init__( self._t_plus = t_plus self._grad_t_plus = grad_t_plus - def proj_kappa_plus(self, dcoll, dd_bdry): + def _normal_kappa_plus(self, dcoll, dd_bdry): # project kappa plus in case of overintegration if isinstance(self._kappa_plus, np.ndarray): # orthotropic materials @@ -266,7 +272,7 @@ def proj_kappa_plus(self, dcoll, dd_bdry): return np.dot(normal, kappa_plus*normal) return project_from_base(dcoll, dd_bdry, self._kappa_plus) - def proj_kappa_minus(self, dcoll, dd_bdry, kappa): + def _normal_kappa_minus(self, dcoll, dd_bdry, kappa): # state minus is already in the quadrature domain if isinstance(kappa, np.ndarray): # orthotropic materials @@ -276,8 +282,8 @@ def proj_kappa_minus(self, dcoll, dd_bdry, kappa): return kappa def kappa_bc(self, dcoll, dd_bdry, kappa_minus): - return harmonic_mean(self.proj_kappa_minus(dcoll, dd_bdry, kappa_minus), - self.proj_kappa_plus(dcoll, dd_bdry)) + return harmonic_mean(kappa_minus, + project_from_base(dcoll, dd_bdry, self._kappa_plus)) def temperature_plus(self, dcoll, dd_bdry): return project_from_base(dcoll, dd_bdry, self._t_plus) @@ -285,9 +291,9 @@ def temperature_plus(self, dcoll, dd_bdry): def temperature_bc(self, dcoll, dd_bdry, kappa_minus, t_minus): t_plus = project_from_base(dcoll, dd_bdry, self._t_plus) actx = t_minus.array_context - kappa_plus = self.proj_kappa_plus(dcoll, dd_bdry) - kappa_minus = self.proj_kappa_minus(dcoll, dd_bdry, - kappa_minus + t_minus*0.0) + kappa_plus = self._normal_kappa_plus(dcoll, dd_bdry) + kappa_minus = self._normal_kappa_minus(dcoll, dd_bdry, + kappa_minus + t_minus*0.0) kappa_sum = actx.np.where( actx.np.greater(kappa_minus + kappa_plus, 0*kappa_minus), kappa_minus + kappa_plus, diff --git a/test/test_diffusion.py b/test/test_diffusion.py index f43bf32f6..08e4d7e10 100644 --- a/test/test_diffusion.py +++ b/test/test_diffusion.py @@ -731,8 +731,7 @@ def test_orthotropic_diffusion(actx_factory): dcoll = create_discretization_collection(actx, mesh, order) nodes = actx.thaw(dcoll.nodes()) - zeros = actx.np.zeros_like(nodes[0]) - kappa = make_obj_array([kappa0 + zeros, kappa1 + zeros]) + kappa = make_obj_array([kappa0, kappa1]) u = 30.0*nodes[0] + 60.0*nodes[1] # exercise Neumann BC @@ -753,8 +752,8 @@ def test_orthotropic_diffusion(actx_factory): assert err_grad_y < 1.e-9 diff_flux = diffusion_flux(kappa, grad_u) - flux_x = -(kappa0 + zeros)*grad_u[0] - flux_y = -(kappa1 + zeros)*grad_u[1] + flux_x = -kappa0*grad_u[0] + flux_y = -kappa1*grad_u[1] err_flux_x = actx.to_numpy(op.norm(dcoll, diff_flux[0] - flux_x, np.inf)) err_flux_y = actx.to_numpy(op.norm(dcoll, diff_flux[1] - flux_y, np.inf)) assert err_flux_x < 1.e-9 @@ -785,8 +784,8 @@ def make_dirichlet_bc(btag): assert err_grad_y < 1.e-9 diff_flux = diffusion_flux(kappa, grad_u) - flux_x = -(kappa0 + zeros)*grad_u[0] - flux_y = -(kappa1 + zeros)*grad_u[1] + flux_x = -kappa0*grad_u[0] + flux_y = -kappa1*grad_u[1] err_flux_x = actx.to_numpy(op.norm(dcoll, diff_flux[0] - flux_x, np.inf)) err_flux_y = actx.to_numpy(op.norm(dcoll, diff_flux[1] - flux_y, np.inf)) assert err_flux_x < 1.e-9 diff --git a/test/test_multiphysics.py b/test/test_multiphysics.py index b00fe4e37..f980389d0 100644 --- a/test/test_multiphysics.py +++ b/test/test_multiphysics.py @@ -224,13 +224,7 @@ def test_thermally_coupled_fluid_wall( # Made-up wall material wall_density = 10*fluid_density wall_heat_capacity = fluid_heat_capacity - if orthotropic_kappa: - wall_kappa = make_obj_array([10*fluid_kappa, 10*fluid_kappa]) - _wall_kappa = wall_kappa[0] # dummy variable to simplify the test - else: - wall_kappa = 10*fluid_kappa - _wall_kappa = wall_kappa # dummy variable to simplify the test - + wall_kappa = 10*fluid_kappa base_wall_temp = 600 fluid_boundaries = { @@ -248,20 +242,20 @@ def test_thermally_coupled_fluid_wall( } interface_temp = ( - (fluid_kappa * base_fluid_temp + _wall_kappa * base_wall_temp) - / (fluid_kappa + _wall_kappa)) + (fluid_kappa * base_fluid_temp + wall_kappa * base_wall_temp) + / (fluid_kappa + wall_kappa)) interface_flux = ( - -fluid_kappa * _wall_kappa / (fluid_kappa + _wall_kappa) + -fluid_kappa * wall_kappa / (fluid_kappa + wall_kappa) * (base_fluid_temp - base_wall_temp)) fluid_alpha = fluid_kappa/(fluid_density * fluid_heat_capacity) - wall_alpha = _wall_kappa/(wall_density * wall_heat_capacity) + wall_alpha = wall_kappa/(wall_density * wall_heat_capacity) def steady_func(kappa, x, t): return interface_temp - interface_flux/kappa * x[1] fluid_steady_func = partial(steady_func, fluid_kappa) - wall_steady_func = partial(steady_func, _wall_kappa) + wall_steady_func = partial(steady_func, wall_kappa) def perturb_func(alpha, x, t): w = 1.5 * np.pi @@ -556,10 +550,98 @@ def test_thermally_coupled_fluid_wall_with_radiation( fluid_rhs = fluid_rhs.energy solid_rhs = wall_energy_rhs/(wall_cp*wall_rho) + # FIXME check convergence instead? assert actx.to_numpy(op.norm(dcoll, fluid_rhs, np.inf)) < 1e-4 assert actx.to_numpy(op.norm(dcoll, solid_rhs, np.inf)) < 1e-4 +@pytest.mark.parametrize("use_overintegration", [False, True]) +@pytest.mark.parametrize("use_noslip", [False, True]) +@pytest.mark.parametrize("use_radiation", [False, True]) +def test_orthotropic_flux( + actx_factory, use_overintegration, use_radiation, use_noslip, + visualize=False): + """Check the RHS shape for orthotropic kappa cases.""" + actx = actx_factory() + + dim = 2 + order = 3 + nelems = 48 + + global_mesh = get_box_mesh(dim, -2, 1, nelems) + + mgrp, = global_mesh.groups + x = global_mesh.vertices[0, mgrp.vertex_indices] + x_elem_avg = np.sum(x, axis=1)/x.shape[0] + volume_to_elements = { + "Fluid": np.where(x_elem_avg > 0)[0], + "Solid": np.where(x_elem_avg <= 0)[0]} + + from meshmode.mesh.processing import partition_mesh + volume_meshes = partition_mesh(global_mesh, volume_to_elements) + + dcoll = create_discretization_collection( + actx, volume_meshes, order=order, quadrature_order=2*order+1) + + quadrature_tag = DISCR_TAG_QUAD if use_overintegration else None + + dd_vol_fluid = DOFDesc(VolumeDomainTag("Fluid"), DISCR_TAG_BASE) + dd_vol_solid = DOFDesc(VolumeDomainTag("Solid"), DISCR_TAG_BASE) + + fluid_nodes = actx.thaw(dcoll.nodes(dd=dd_vol_fluid)) + solid_nodes = actx.thaw(dcoll.nodes(dd=dd_vol_solid)) + + eos = IdealSingleGas() + transport = SimpleTransport(viscosity=0.0, thermal_conductivity=10.0) + gas_model = GasModel(eos=eos, transport=transport) + + # Fluid cv + + mom_x = 0.0*fluid_nodes[0] + mom_y = 0.0*fluid_nodes[0] + momentum = make_obj_array([mom_x, mom_y]) + + temperature = 2.0 + 0.0*fluid_nodes[0] + mass = 1.0 + 0.0*fluid_nodes[0] + energy = mass*eos.get_internal_energy(temperature) + + fluid_cv = make_conserved(dim=dim, mass=mass, energy=energy, + momentum=momentum) + + fluid_state = make_fluid_state(cv=fluid_cv, gas_model=gas_model) + + # Made-up wall material + wall_kappa = make_obj_array([1.0, 1.0]) + wall_emissivity = 1.0 + wall_temperature = 1.0 + 0.0*solid_nodes[0] + + fluid_boundaries = { + dd_vol_fluid.trace("-2").domain_tag: AdiabaticNoslipWallBoundary(), + dd_vol_fluid.trace("+2").domain_tag: AdiabaticNoslipWallBoundary(), + dd_vol_fluid.trace("+1").domain_tag: + IsothermalWallBoundary(wall_temperature=2.0), + } + + solid_boundaries = { + dd_vol_solid.trace("-2").domain_tag: NeumannDiffusionBoundary(0.), + dd_vol_solid.trace("+2").domain_tag: NeumannDiffusionBoundary(0.), + dd_vol_solid.trace("-1").domain_tag: DirichletDiffusionBoundary(1.0), + } + + fluid_rhs, wall_energy_rhs = coupled_ns_heat_operator( + dcoll, gas_model, dd_vol_fluid, dd_vol_solid, fluid_boundaries, + solid_boundaries, fluid_state, wall_kappa, wall_temperature, + time=0.0, quadrature_tag=quadrature_tag, interface_noslip=use_noslip, + interface_radiation=use_radiation, sigma=2.0, + ambient_temperature=0.0, wall_emissivity=wall_emissivity, + inviscid_terms_on=False + ) + + from meshmode.dof_array import DOFArray + assert isinstance(fluid_rhs.energy, DOFArray) + assert isinstance(wall_energy_rhs, DOFArray) + + if __name__ == "__main__": import sys if len(sys.argv) > 1: