|
| 1 | +# Import required libs |
| 2 | +from fenics import Constant, Function, AutoSubDomain, VectorFunctionSpace, interpolate, \ |
| 3 | + TrialFunction, TestFunction, Point, Expression, DirichletBC, nabla_grad, \ |
| 4 | + Identity, inner, dx, ds, sym, grad, lhs, rhs, dot, File, solve, assemble_system |
| 5 | +from mshr import Cylinder, generate_mesh |
| 6 | +from ufl import nabla_div |
| 7 | +import numpy as np |
| 8 | +from fenicsprecice import Adapter |
| 9 | +import math |
| 10 | + |
| 11 | + |
| 12 | +# define the two kinds of boundary: clamped and coupling Neumann Boundary |
| 13 | +def clamped_boundary(x, on_boundary): |
| 14 | + """ |
| 15 | + Filter nodes at both ends of tube as they are fixed |
| 16 | + """ |
| 17 | + tol = 1E-14 |
| 18 | + return on_boundary and (((abs(x[2]) - 0.0) < tol) or ((L - abs(x[2])) < tol)) |
| 19 | + |
| 20 | + |
| 21 | +def neumann_boundary(x, on_boundary): |
| 22 | + """ |
| 23 | + Filter nodes which lie on the inner surface of the tube and excluding end nodes |
| 24 | + """ |
| 25 | + tol = 1E-14 |
| 26 | + return on_boundary and ((math.sqrt(x[0]**2 + x[1]**2) - R) < tol) and ((L - x[2]) > tol) and ((x[2] - 0.0) > tol) |
| 27 | + |
| 28 | + |
| 29 | +# Geometry and material properties |
| 30 | +dim = 3 # number of dimensions |
| 31 | +R = 0.005 |
| 32 | +L = 0.05 |
| 33 | +rho = 3000 |
| 34 | +E = 4000000 |
| 35 | +nu = 0.3 |
| 36 | + |
| 37 | +mu = Constant(E / (2.0 * (1.0 + nu))) |
| 38 | + |
| 39 | +lambda_ = Constant(E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))) |
| 40 | + |
| 41 | +# create Mesh |
| 42 | +outer_tube = Cylinder(Point(0, 0, L), Point(0, 0, 0), R + 0.001, R + 0.001) |
| 43 | +inner_tube = Cylinder(Point(0, 0, L), Point(0, 0, 0), R, R) |
| 44 | +mesh = generate_mesh(outer_tube - inner_tube, 20) |
| 45 | + |
| 46 | +# create Function Space |
| 47 | +V = VectorFunctionSpace(mesh, 'P', 2) |
| 48 | + |
| 49 | +# Trial and Test Functions |
| 50 | +du = TrialFunction(V) |
| 51 | +v = TestFunction(V) |
| 52 | + |
| 53 | +u_np1 = Function(V) |
| 54 | +saved_u_old = Function(V) |
| 55 | +u_delta = Function(V) |
| 56 | + |
| 57 | +# function known from previous timestep |
| 58 | +u_n = Function(V) |
| 59 | +v_n = Function(V) |
| 60 | +a_n = Function(V) |
| 61 | + |
| 62 | +coupling_boundary = AutoSubDomain(neumann_boundary) |
| 63 | +fixed_boundary = AutoSubDomain(clamped_boundary) |
| 64 | + |
| 65 | +precice = Adapter(adapter_config_filename="precice-adapter-config-fsi-s.json") |
| 66 | + |
| 67 | +# Initialize the coupling interface |
| 68 | +precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=V, fixed_boundary=fixed_boundary) |
| 69 | + |
| 70 | +fenics_dt = precice_dt # if fenics_dt == precice_dt, no subcycling is applied |
| 71 | +dt = Constant(np.min([precice_dt, fenics_dt])) |
| 72 | + |
| 73 | +# clamp the tube on both sides |
| 74 | +bc = DirichletBC(V, Constant((0, 0, 0)), fixed_boundary) |
| 75 | + |
| 76 | +# alpha method parameters |
| 77 | +alpha_m = Constant(0) |
| 78 | +alpha_f = Constant(0) |
| 79 | +gamma = Constant(0.5 + alpha_f - alpha_m) |
| 80 | +beta = Constant((gamma + 0.5) ** 2 / 4.) |
| 81 | + |
| 82 | + |
| 83 | +# Define strain |
| 84 | +def epsilon(u): |
| 85 | + return 0.5 * (nabla_grad(u) + nabla_grad(u).T) |
| 86 | + |
| 87 | + |
| 88 | +# Define Stress tensor |
| 89 | +def sigma(u): |
| 90 | + return lambda_ * nabla_div(u) * Identity(dim) + 2 * mu * epsilon(u) |
| 91 | + |
| 92 | + |
| 93 | +# Define Mass form |
| 94 | +def m(u, v): |
| 95 | + return rho * inner(u, v) * dx |
| 96 | + |
| 97 | + |
| 98 | +# Elastic stiffness form |
| 99 | +def k(u, v): |
| 100 | + return inner(sigma(u), sym(grad(v))) * dx |
| 101 | + |
| 102 | + |
| 103 | +# External Work |
| 104 | +def Wext(u_): |
| 105 | + return dot(u_, p) * ds |
| 106 | + |
| 107 | + |
| 108 | +# Functions for updating system state |
| 109 | + |
| 110 | +# Update acceleration |
| 111 | +def update_a(u, u_old, v_old, a_old, ufl=True): |
| 112 | + if ufl: |
| 113 | + dt_ = dt |
| 114 | + beta_ = beta |
| 115 | + else: |
| 116 | + dt_ = float(dt) |
| 117 | + beta_ = float(beta) |
| 118 | + |
| 119 | + return ((u - u_old - dt_ * v_old) / beta / dt_ ** 2 |
| 120 | + - (1 - 2 * beta_) / 2 / beta_ * a_old) |
| 121 | + |
| 122 | + |
| 123 | +# Update velocity |
| 124 | +def update_v(a, u_old, v_old, a_old, ufl=True): |
| 125 | + if ufl: |
| 126 | + dt_ = dt |
| 127 | + gamma_ = gamma |
| 128 | + else: |
| 129 | + dt_ = float(dt) |
| 130 | + gamma_ = float(gamma) |
| 131 | + |
| 132 | + return v_old + dt_ * ((1 - gamma_) * a_old + gamma_ * a) |
| 133 | + |
| 134 | + |
| 135 | +def update_fields(u, u_old, v_old, a_old): |
| 136 | + """Update all fields at the end of a timestep.""" |
| 137 | + |
| 138 | + u_vec, u0_vec = u.vector(), u_old.vector() |
| 139 | + v0_vec, a0_vec = v_old.vector(), a_old.vector() |
| 140 | + |
| 141 | + # call update functions |
| 142 | + a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False) |
| 143 | + v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False) |
| 144 | + |
| 145 | + # assign u->u_old |
| 146 | + v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec |
| 147 | + u_old.vector()[:] = u.vector() |
| 148 | + |
| 149 | + |
| 150 | +def avg(x_old, x_new, alpha): |
| 151 | + return alpha * x_old + (1 - alpha) * x_new |
| 152 | + |
| 153 | + |
| 154 | +a_np1 = update_a(du, u_n, v_n, a_n, ufl=True) |
| 155 | +v_np1 = update_v(a_np1, u_n, v_n, a_n, ufl=True) |
| 156 | + |
| 157 | +res = m(avg(a_n, a_np1, alpha_m), v) + k(avg(u_n, du, alpha_f), v) |
| 158 | + |
| 159 | +a_form = lhs(res) |
| 160 | +L_form = rhs(res) |
| 161 | + |
| 162 | +# parameters for Time-Stepping |
| 163 | +t = 0.0 |
| 164 | +n = 0 |
| 165 | +E_ext = 0 |
| 166 | + |
| 167 | +displacement_out = File("output/u_fsi.pvd") |
| 168 | + |
| 169 | +u_n.rename("Displacement", "") |
| 170 | +u_np1.rename("Displacement", "") |
| 171 | +displacement_out << u_n |
| 172 | + |
| 173 | +while precice.is_coupling_ongoing(): |
| 174 | + |
| 175 | + if precice.is_action_required(precice.action_write_iteration_checkpoint()): # write checkpoint |
| 176 | + precice.store_checkpoint(u_n, t, n) |
| 177 | + |
| 178 | + # read data from preCICE and get a new coupling expression |
| 179 | + read_data = precice.read_data() |
| 180 | + |
| 181 | + # Update the point sources on the coupling boundary with the new read data |
| 182 | + forces_x, forces_y, forces_z = precice.get_point_sources(read_data) |
| 183 | + |
| 184 | + A, b = assemble_system(a_form, L_form, bc) |
| 185 | + |
| 186 | + b_forces = b.copy() # b is the same for every iteration, only forces change |
| 187 | + |
| 188 | + for ps in forces_x: |
| 189 | + ps.apply(b_forces) |
| 190 | + for ps in forces_y: |
| 191 | + ps.apply(b_forces) |
| 192 | + for ps in forces_z: |
| 193 | + ps.apply(b_forces) |
| 194 | + |
| 195 | + assert (b is not b_forces) |
| 196 | + solve(A, u_np1.vector(), b_forces) |
| 197 | + |
| 198 | + dt = Constant(np.min([precice_dt, fenics_dt])) |
| 199 | + |
| 200 | + # Write relative displacements to preCICE |
| 201 | + u_delta.vector()[:] = u_np1.vector()[:] - u_n.vector()[:] |
| 202 | + precice.write_data(u_delta) |
| 203 | + |
| 204 | + # Call to advance coupling, also returns the optimum time step value |
| 205 | + precice_dt = precice.advance(dt(0)) |
| 206 | + |
| 207 | + # Either revert to old step if timestep has not converged or move to next timestep |
| 208 | + if precice.is_action_required(precice.action_read_iteration_checkpoint()): # roll back to checkpoint |
| 209 | + u_cp, t_cp, n_cp = precice.retrieve_checkpoint() |
| 210 | + u_n.assign(u_cp) |
| 211 | + t = t_cp |
| 212 | + n = n_cp |
| 213 | + else: |
| 214 | + u_n.assign(u_np1) |
| 215 | + t += float(dt) |
| 216 | + n += 1 |
| 217 | + |
| 218 | + if precice.is_time_window_complete(): |
| 219 | + update_fields(u_np1, saved_u_old, v_n, a_n) |
| 220 | + if n % 10 == 0: |
| 221 | + displacement_out << (u_n, t) |
| 222 | + |
| 223 | +# Plot tip displacement evolution |
| 224 | +displacement_out << u_n |
| 225 | + |
| 226 | +precice.finalize() |
0 commit comments