From 3b2cecfe22bb9a7c1af64ca50d0bc60f16c0eee2 Mon Sep 17 00:00:00 2001 From: Martin Bauer <martin.bauer@fau.de> Date: Mon, 3 Sep 2018 12:59:37 +0200 Subject: [PATCH] lbm phasefield improvements - 1D benchmark scenario, trying different correction functions - 2D neumann angle evaluation based on circle fitting --- phasefield/analytical.py | 6 ++ phasefield/contact_angle_circle_fitting.py | 120 +++++++++++++++++++++ phasefield/phasefieldstep.py | 16 ++- phasefield/scenarios.py | 1 + 4 files changed, 141 insertions(+), 2 deletions(-) create mode 100644 phasefield/contact_angle_circle_fitting.py diff --git a/phasefield/analytical.py b/phasefield/analytical.py index 2a6ef6e4..2a3854e2 100644 --- a/phasefield/analytical.py +++ b/phasefield/analytical.py @@ -78,6 +78,12 @@ def free_energy_functional_n_phases_penalty_term(order_parameters, interface_wid return bulk + interface + penalty_term_factor * bulk_penalty_term +def n_phases_correction_function(c, beta, power=2): + return sp.Piecewise((-beta * c ** power, c < 0), + (-beta * (1 - c) ** 2, c > 1), + (c ** 2 * (1 - c) ** power, True)) + + def free_energy_functional_n_phases(num_phases=None, surface_tensions=symmetric_symbolic_surface_tension, interface_width=interface_width_symbol, order_parameters=None, include_bulk=True, include_interface=True, symbolic_lambda=False, diff --git a/phasefield/contact_angle_circle_fitting.py b/phasefield/contact_angle_circle_fitting.py new file mode 100644 index 00000000..fa70e220 --- /dev/null +++ b/phasefield/contact_angle_circle_fitting.py @@ -0,0 +1,120 @@ +from collections import namedtuple + +import numpy as np + + +def circle_intersections(midpoint0, midpoint1, radius0, radius1): + """Returns intersection points of two circles.""" + midpoint_distance = np.linalg.norm(midpoint0 - midpoint1) + x0, y0 = midpoint0 + x1, y1 = midpoint1 + + k1 = (y0 - y1) / (x0 - x1) # slope for the line linking two centers + b1 = y1 - k1 * x1 # line equation of the line linking two centers + + k2 = -1 / k1 # slope for the line linking two cross points + b2 = (radius0 ** 2 - radius1 ** 2 - x0 ** 2 + x1 ** 2 - y0 ** 2 + y1 ** 2) / (2 * (y1 - y0)) + + if midpoint_distance > (radius0 + radius1): # no intersection + return [] + + if np.abs(radius1 - radius0) == midpoint_distance or midpoint_distance == radius0 + radius1: + xx = -(b1 - b2) / (k1 - k2) + yy = -(-b2 * k1 + b1 * k2) / (k1 - k2) + return [(xx, yy)] + elif np.abs(radius1 - radius0) < midpoint_distance < radius0 + radius1: + xx1 = (-b2 * k2 + x1 + k2 * y1 - np.sqrt(-b2 ** 2 + radius1 ** 2 + k2 ** 2 * radius1 ** 2 - 2 * b2 * + k2 * x1 - k2 ** 2 * x1 ** 2 + 2 * b2 * y1 + 2 * k2 * x1 * y1 - y1 ** 2)) / (1 + k2 ** 2) + yy1 = k2 * xx1 + b2 + + xx2 = (-b2 * k2 + x1 + k2 * y1 + np.sqrt(-b2 ** 2 + radius1 ** 2 + k2 ** 2 * radius1 ** 2 - 2 * b2 * + k2 * x1 - k2 ** 2 * x1 ** 2 + 2 * b2 * y1 + 2 * k2 * x1 * y1 - y1 ** 2)) / (1 + k2 ** 2) + yy2 = k2 * xx2 + b2 + + return [(xx1, yy1), (xx2, yy2)] + else: + return [] + + +def interface_region(concentration_arr, phase0, phase1, area=3): + import scipy.ndimage as sc_image + + area_phase0 = sc_image.morphology.binary_dilation(concentration_arr[..., phase0] > 0.5, iterations=area) + area_phase1 = sc_image.morphology.binary_dilation(concentration_arr[..., phase1] > 0.5, iterations=area) + return np.logical_and(area_phase0, area_phase1) + + +def contour_line(concentration_arr, phase0, phase1, cutoff=0.05): + from skimage import measure + concentration_arr = concentration_arr.copy() + + mask = interface_region(concentration_arr, phase0, phase1) + concentration_arr = concentration_arr[..., phase0] + concentration_arr[np.logical_not(mask)] = np.nan + contours = measure.find_contours(concentration_arr, 0.5) + + contours = [c for c in contours if len(c) > 8] + + if len(contours) != 1: + raise ValueError("Multiple or non contour lines (%d) found" % (len(contours),)) + + contour = contours[0] + absolute_cutoff = int(len(contour) * cutoff) + 1 + return contour[absolute_cutoff:-absolute_cutoff] + + +def fit_circle(points): + # see https://scipy-cookbook.readthedocs.io/items/Least_Squares_Circle.html + from scipy import optimize + + def point_distances(xc, yc): + """ calculate the distance of each 2D points from the center (xc, yc) """ + return np.sqrt((points[:, 0]-xc)**2 + (points[:, 1]-yc)**2) + + def f(c): + """ calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc) """ + ri = point_distances(*c) + return ri - ri.mean() + + center_estimate = np.mean(points, axis=0) + center, _ = optimize.leastsq(f, center_estimate) + radius = point_distances(*center).mean() + + circle = namedtuple('Circle', ['center', 'radius']) + return circle(center, radius) + + +def neumann_angles_from_surface_tensions(surface_tensions): + g = surface_tensions + angles = [ + np.arccos(-(g(0, 1) * g(0, 1) + g(0, 2) * g(0, 2) - g(1, 2) * g(1, 2)) / 2 / g(0, 1) / g(0, 2)), + np.arccos(-(g(0, 1) * g(0, 1) + g(1, 2) * g(1, 2) - g(0, 2) * g(0, 2)) / 2 / g(0, 1) / g(1, 2)), + np.arccos(-(g(1, 2) * g(1, 2) + g(0, 2) * g(0, 2) - g(0, 1) * g(0, 1)) / 2 / g(1, 2) / g(0, 2)), + ] + return [np.rad2deg(a) for a in angles] + + +def liquid_lens_neumann_angles(concentration, drop_phase_idx=2, enclosing_phase1=0, enclosing_phase2=1): + circles = [fit_circle(contour_line(concentration, enclosing_phase1, drop_phase_idx)), + fit_circle(contour_line(concentration, enclosing_phase2, drop_phase_idx))] + + intersections = circle_intersections(circles[0].center, circles[1].center, circles[0].radius, circles[1].radius) + + y1, y2 = (i[1] for i in intersections) + assert np.abs(y1 - y2) < 0.1, "Liquid lens is not aligned in y direction (rotated?)" + y_horizontal = (y1 + y2) / 2 + + xs = np.sqrt(circles[0].radius ** 2 - (y_horizontal - circles[0].center[1]) ** 2) + angle = np.rad2deg(np.arctan(-xs / np.sqrt(circles[0].radius ** 2 - xs ** 2))) % 180 + angle = 180 - angle + neumann3_upper = angle + neumann1 = 180 - neumann3_upper + + xs = np.sqrt(circles[1].radius ** 2 - (y_horizontal - circles[1].center[1]) ** 2) + angle = np.rad2deg(np.arctan(-xs / np.sqrt(circles[1].radius ** 2 - xs ** 2))) % 180 + angle = 180 - angle + neumann3_lower = angle + neumann2 = 180 - neumann3_lower + + neumann3 = neumann3_upper + neumann3_lower + return neumann1, neumann2, neumann3 diff --git a/phasefield/phasefieldstep.py b/phasefield/phasefieldstep.py index a5f48169..45a0fba8 100644 --- a/phasefield/phasefieldstep.py +++ b/phasefield/phasefieldstep.py @@ -169,7 +169,6 @@ class PhaseFieldStep: if op == density_order_parameter: continue - print("CH Gamma", cahn_hilliard_gammas[i]) ch_method = cahn_hilliard_lb_method(self.hydro_lbm_step.method.stencil, self.mu_field(i), relaxation_rate=cahn_hilliard_relaxation_rates[i], gamma=cahn_hilliard_gammas[i]) @@ -295,8 +294,21 @@ class PhaseFieldStep: return self._get_slice(self.phi_field_name, slice_obj) def concentration_slice(self, slice_obj=None): + if slice_obj is not None and len(slice_obj) > self.data_handling.dim: + assert len(slice_obj) - 1 == self.data_handling.dim + last_index = slice_obj[-1] + slice_obj = slice_obj[:-1] + else: + last_index = None + phi = self.phi_slice(slice_obj) - return phi if self.order_parameters_to_concentrations is None else self.order_parameters_to_concentrations(phi) + if self.order_parameters_to_concentrations is None: + return phi + else: + result = self.order_parameters_to_concentrations(phi) + if last_index is not None: + result = result[..., last_index] + return result def mu_slice(self, slice_obj=None): return self._get_slice(self.mu_field_name, slice_obj) diff --git a/phasefield/scenarios.py b/phasefield/scenarios.py index 851446d1..807970ea 100644 --- a/phasefield/scenarios.py +++ b/phasefield/scenarios.py @@ -36,6 +36,7 @@ def create_three_phase_model(alpha=1, kappa=(0.015, 0.015, 0.015), include_rho=T return PhaseFieldStep(free_energy, order_parameters, density_order_parameter=None, transformation_matrix=transformation_matrix, order_parameters_to_concentrations=order_parameters_to_concentrations, + cahn_hilliard_gammas=[1, 1, 1 / 3], **kwargs) -- GitLab