Add load over a surfaces in a rotatory geometry

dolfinx
, , , ,
1

Hi, I’m new to FenicsX, so excuse me if this is too basic. I’m coming from this tutorial trying to simulate the elastic response of a wind turbine rotor, so I’m dealing with a transient elastic problem.

The problem involves a two-blade rotor, where I want to consider gravity, centrifugal force, and wind speed load. The rotor rotates around the X-axis. Currently, I’m using a fixed vector as a wind load, but eventually, this load will be derived from real data or CFD computations. For now, it’s a single vector.

Right now, I have two questions:

  1. How do I apply the wind load only over the rotor surface? (The mesh and geometry are created with Gmsh and read from the .msh file.), At this point I discover how to apply it over all the surface, but I want to apply over the selected surface marked with the tag GMSH_TAG_SURFACE. Searching on the forum I found several answers about using subdomains like this:
ds = fe.Measure('ds', domain=mesh, subdomain_data=subdomains)

but I can’t find how to create the subdomain in my example, with the FenicsX

  1. I’m not sure if I’m heading in the right direction by handling the rotation through applying geometry transformations in the solve loop. I’m uncertain if this breaks the continuity of the system. I think that I have to apply the rotation also to the displacment field, Im really lost here

Here is my code at the moment:
Mesh can be downloaded here: Very Coarse Mesh File

I have modified my example code to give a working example and fix some errors…

import numpy as np

from mpi4py import MPI
import ufl
import dolfinx.fem.petsc
from dolfinx import fem, io, default_scalar_type

from scipy.spatial.transform import Rotation

comm = MPI.COMM_WORLD

# ------------- #
# Problem Setup #
# ------------- #
# material properties
E = 210e3
nu = 0.3
rho = 2e-3

# fixed rotation speed
rpm = -72  # rpm
omega = (rpm / 60) * 2 * ufl.pi

# loads
gravity = (0, 9.81, 0)
wind_load = (10, 0, 0)

# Timing
startTime = 0.0
endTime = 30.0
deltaT = 0.001
writeInterval = 1 # timestep

# Newmark parameters
# which ensures unconditional stability,
# energy conservation and second-order accuracy
beta_ = 0.25
gamma_ = 0.5

# ------- #
# MESHING #
# ------- #
GMSH_TAG_FIX = 1003
GMSH_TAG_SURFACE = 1004
mesh_file = "NREL_Phase_VI_rotor.msh"  # attached file url
domain, cell_markers, facet_markers = io.gmshio.read_from_msh(mesh_file, comm, gdim=3)

# -------------- #
# Function Space #
# -------------- #
dim = domain.geometry.dim
degree = 2
V = fem.functionspace(domain, ("Lagrange", degree, (dim,)))
u = fem.Function(V, name="Displacement")

# Here how to select only the GMSH_TAG_SURFACE as subdomain
# for surface integration ?
ds = ufl.Measure("ds", domain=domain)

# ------------------- #
# Boundary conditions #
# ------------------- #
# Wind Load
# rotor_facets = facet_markers.find(GMSH_TAG_SURFACE)
# rotor_surf = fem.locate_dofs_topological(V, 2, rotor_facets)
# f_w = fem.Constant(domain, default_scalar_type(wind_load))  # <== how apply only to the rotor_facets?

# Gravitational load
f_g = fem.Constant(domain, default_scalar_type(gravity))

# Centrifugal acceleration (-rho*omega^2*x_i) load
# x = ufl.SpatialCoordinate(domain)
# f_c = ufl.as_vector((0, -rho * omega**2 * x[1], -rho * omega**2 * x[2]))

# Rotation angle
theta = fem.Constant(domain, 0.0)

# Fixed BC
fixed_facets = facet_markers.find(GMSH_TAG_FIX)
fixed_surface = fem.locate_dofs_topological(V, 2, fixed_facets)
fixed_bc = fem.dirichletbc(np.array([0.0, 0.0, 0.0]), fixed_surface, V)


# ------------------------ #
# linear elastic equations #
# ------------------------ #
E = fem.Constant(domain, E)
nu = fem.Constant(domain, nu)
rho = fem.Constant(domain, rho)

lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2 / (1 + nu)


def epsilon(v):
    return ufl.sym(ufl.grad(v))


def sigma(v):
    return lmbda * ufl.tr(epsilon(v)) * ufl.Identity(dim) + 2 * mu * epsilon(v)


# ------------------- #
# Time discretization #
# ------------------- #
# prev time step
u_old = fem.Function(V)
v_old = fem.Function(V)
a_old = fem.Function(V)

# current time step
a_new = fem.Function(V)
v_new = fem.Function(V)

beta = fem.Constant(domain, beta_)
gamma = fem.Constant(domain, gamma_)
dt = fem.Constant(domain, 0.0)

a = 1 / beta / dt**2 * (u - u_old - dt * v_old) + a_old * (1 - 1 / 2 / beta)
a_expr = fem.Expression(a, V.element.interpolation_points())

v = v_old + dt * ((1 - gamma) * a_old + gamma * a)
v_expr = fem.Expression(v, V.element.interpolation_points())

# ------------------------------- #
# mass, a stiffness and a damping #
# ------------------------------- #
u_ = ufl.TestFunction(V)
du = ufl.TrialFunction(V)

eta_M = fem.Constant(domain, 1e-4)
eta_K = fem.Constant(domain, 1e-4)


def mass(u, u_):
    return rho * ufl.dot(u, u_) * ufl.dx


def stiffness(u, u_):
    return ufl.inner(sigma(u), epsilon(u_)) * ufl.dx


def damping(u, u_):
    return eta_M * mass(u, u_) + eta_K * stiffness(u, u_)


L = ufl.dot(f_g, u_) * ufl.dx
# L += ufl.dot(f_c, u_) * ufl.dx # include centrigual force
# L += ufl.dot(f_w, u_) * ufl.ds # include wind loads over surface
Residual = mass(a, u_) + damping(v, u_) + stiffness(u, u_) - L

Residual_du = ufl.replace(Residual, {u: du})
a_form = ufl.lhs(Residual_du)
L_form = ufl.rhs(Residual_du)

bcs = [fixed_bc]

problem = fem.petsc.LinearProblem(
    a_form, L_form, u=u, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)

# ---------- #
# Solve loop #
# ---------- #
time_steps = int(endTime / deltaT)
t = startTime
i = 0
current_theta = 0
vtk = io.VTKFile(domain.comm, "test_results/elastodynamics.pvd", "w")
running = True


while running:
    if t > 0:
        # start rotation arown X axis, this is correct?
        theta.value = omega * deltaT
        rotation = Rotation.from_euler("x", theta.value, degrees=False)
        domain.geometry.x[:, :] = rotation.apply(domain.geometry.x)
        current_theta += theta.value

    if i % writeInterval == 0:
        vtk.write_function(u, t)

    dt.value = deltaT
    t += deltaT

    # do I have to re-apply the loads here due to the domain rotatio?
    f_g.value = np.array(gravity)
    # f_w.value = np.array(wind_load)

    problem.solve()

    u.x.scatter_forward()  # updates ghost values for parallel computations

    # compute new velocity v_n+1
    v_new.interpolate(v_expr)

    # compute new acceleration a_n+1
    a_new.interpolate(a_expr)

    # update u_n with u_n+1
    u.vector.copy(u_old.vector)

    # update v_n with v_n+1
    v_new.vector.copy(v_old.vector)

    # update a_n with a_n+1
    a_new.vector.copy(a_old.vector)
    if comm.rank == 0:
        print(f"Time: {t:.0e}s, Current Angle: {np.rad2deg(current_theta)%360:.2f}º")

    running = t <= endTime
    i += 1

vtk.close()
2

From your call to gmshio

domain, cell_markers, facet_markers = io.gmshio.read_from_msh(mesh_file, comm, gdim=3)

the right way to construct the ds measure is

ds = ufl.Measure('ds', subdomain_data=facet_markers)

Please elaborate why you think that is happening. If by system you mean “mesh”, my understanding from the code is that you are rotating the entire geometry, so there shouldn’t be any tear in the mesh after rotation.

1 Like
3

Dear @francesco-ballarin thanks for your answer. So If, I want to create a integration surface only for a certain part of the domain, I can do something like this:

ds = ufl.Measure("ds", subdomain_data=facet_markers.find(GMSH_TAG_SURFACE))

And doing this I avoid the integration over the GMSH_TAG_FIX surface, which represent a fixed part of the domain, and it’s not under the wind loads effect, in my case. Is that correct?.

When I’m talking abot breaking the contiuity of the system, is that… for example. As I’m roating the mesh coordinates, after defining the model, and I’m doing it “by hand” I’m not sure if that angular displacement that I’m introdicing with the rotation is taken into account for system solution or if its just solving every new position as a snapshoot of the system in space/time, without taking into account the angular inertial moment. That’s why I’m introducing de centripetal force, but, I’m not sure if that is sufficient.

4

If you want to integrate over a single surface, use

ds = ufl.Measure('ds', subdomain_data=facet_markers)
my_form = v * ds(GMSH_TAG_SURFACE)

As I’m roating the mesh coordinates

I can only comment on the technical detail on how you update the mesh coordinates. From a technical standpoint, updating domain.geometry.x looks correct to me. I am not familiar enough with what you mention in the rest of the sentence to comment any further.

1 Like
5

@francesco-ballarin Thanks a lot for your help. I think that I’m gin a run some simulaion and see the results to see If the physics is rasonable after applying the rotation.