How to apply point forces on corners of a 3D plate?

variational formulation
1

Hi guys,

I am tried to perform an eigenanalysis of a plate where tensile or spring loads applied at the bottom four corners of the plate.

Can anyone tell me how to apply these loads. I read similar topics for guidance, but I cannot find one that applies to my case. Here is the code so far.

# importing libraries for pre-processing

import dolfin as df
from ufl import nabla_div, nabla_grad
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
import json
from mshr import Polygon, generate_mesh
from scipy.optimize import root
from math import cos, cosh
import sys
import argparse
from functools import partial
import pandas as pd
import meshio as mio
from stl import mesh as mstl
np.set_printoptions(precision=5)

## Plate Parameters

# Universal Constants
GRAVITY = 9.81
WATER_DENSITY = 1000
PLATE_LENGTH = 1
WATER_DEPTH = 1

# Plate Data
PLATE_DEPTH = 0.2 * WATER_DEPTH
PLATE_WIDTH = 16 * PLATE_LENGTH
PLATE_THICKNESS = 1e-3 * PLATE_LENGTH
PLATE_DENSITY = 1 * WATER_DENSITY
YOUNGS_MODULUS = 1e6 * (WATER_DENSITY * GRAVITY * WATER_DEPTH**4)
POISSONS_RATIO = 0.3
LAMES_MU = 0.5 * YOUNGS_MODULUS / (1+POISSONS_RATIO)
LAMES_LAMBDA = YOUNGS_MODULUS * POISSONS_RATIO / ((1+POISSONS_RATIO)*(1-2*POISSONS_RATIO))

# Tension load at the bottom four corners of the plate
TENSION_LOADS = np.array([0.05, 0.15, 1e4]) * (WATER_DENSITY * GRAVITY * WATER_DEPTH**4)

# Mode Data
MODE_COUNT = 6

# Meshing Data
PANEL_RES = np.array([20, 60, 3],dtype=int)
PANEL_HIGH_RES = np.array([20, 160, 6],dtype=int)

# Simulation directories
MAIN_DIR = Path(".")
OUTPUT_DIR = MAIN_DIR / "output"
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
PLOTS_DIR = MAIN_DIR / "plots"
PLOTS_DIR.mkdir(parents=True, exist_ok=True)

## Meshing

PLATE_THICKNESS/PANEL_RES[2]

CORNER1 = df.Point(0, 0, -PLATE_THICKNESS-PLATE_DEPTH)

# corner2 is diagonally opposite to the corner1
CORNER2 = df.Point(PLATE_LENGTH, 
                   PLATE_WIDTH, 
                   -PLATE_DEPTH)

# res is the resolution of the mesh
low_res_mesh = df.BoxMesh(CORNER1, CORNER2, PANEL_RES[0], PANEL_RES[1], PANEL_RES[2])

# Saving the mesh
df.File("low_res_mesh.pvd") << low_res_mesh
df.File("low_res_mesh.xml") << low_res_mesh

mesh = df.BoxMesh(CORNER1, CORNER2, 
                          PANEL_HIGH_RES[0], PANEL_HIGH_RES[1], PANEL_HIGH_RES[2])

# Saving the mesh
df.File("hig_res_mesh.pvd") << mesh
df.File("high_res_mesh.xml") << mesh
DIM = mesh.geometric_dimension()

MESH_SIZES = np.array([PLATE_LENGTH/PANEL_HIGH_RES[0], 
                    PLATE_WIDTH/PANEL_HIGH_RES[1], 
                    PLATE_THICKNESS/PANEL_HIGH_RES[2]])
print(f"{MESH_SIZES=}")
mesh







# Saving the mesh
df.File("mesh.pvd") << mesh
df.File("mesh.xml") << mesh

## Strain and Stress Functions

# Strain Function
def strain(disp):
    return df.sym(df.grad(disp))

# Stress Function
def stress(disp):
    dim = disp.geometric_dimension()
    return 2.0*LAMES_MU*strain(disp) + \
           LAMES_LAMBDA*df.tr(strain(disp))*df.Identity(dim)

## Applying Constraint Boundary Conditions at the corners of the plate bottom
Note: Following code also has option to add boundary condition at the corners

V = df.VectorFunctionSpace(mesh, 'Lagrange', degree=1)
u_ = df.TrialFunction(V)
du = df.TestFunction(V)

# https://fenicsproject.org/qa/9049/how-to-create-a-dirichletbc-on-a-single-vertex-or-dof/
def corner_spring_bc(x, on_boundary):
    """Function to apply spring boundary condition at the corners of the plate bottom"""
    bottom = df.near(x[2],-PLATE_THICKNESS,MESH_SIZES[2]) and on_boundary
    
    if bottom:
        left   = df.near(x[0], 0,            MESH_SIZES[0])
        right  = df.near(x[0], PLATE_LENGTH, MESH_SIZES[0])
        front  = df.near(x[1], PLATE_WIDTH,  MESH_SIZES[1])
        back   = df.near(x[1], 0,            MESH_SIZES[1])

        corner1 = left  and front
        corner2 = left  and back
        corner3 = right and front
        corner4 = right and back
    
        return corner1 or corner2 or corner3 or corner4
    
    return bottom

def right_corners(x, on_boundary):
    """Function to apply spring boundary condition at the corners of the plate bottom"""
    bottom = df.near(x[2],-PLATE_THICKNESS,MESH_SIZES[2]) and on_boundary
    
    if bottom:
        right  = df.near(x[0], PLATE_LENGTH, MESH_SIZES[0])
        front  = df.near(x[1], PLATE_WIDTH,  MESH_SIZES[1])
        back   = df.near(x[1], 0,            MESH_SIZES[1])

        corner3 = right and front
        corner4 = right and back
    
        return corner3 or corner4
    
    return bottom

def left_corners(x, on_boundary):
    """Function to apply spring boundary condition at the corners of the plate bottom"""
    bottom = df.near(x[2],-PLATE_THICKNESS,MESH_SIZES[2]) and on_boundary
    
    if bottom:
        left   = df.near(x[0], 0,            MESH_SIZES[0])
        front  = df.near(x[1], PLATE_WIDTH,  MESH_SIZES[1])
        back   = df.near(x[1], 0,            MESH_SIZES[1])

        corner1 = left  and front
        corner2 = left  and back
        
        return corner1 or corner2
    
    return bottom

# bc = df.DirichletBC(V.sub(0), 0, corner_spring_bc, method='pointwise')
bc = None

# Applying Tensile Loads

corner1_load = df.PointSource(V.sub(0), df.Point(0,           0          ,0),TENSION_LOADS[0])
corner2_load = df.PointSource(V.sub(0), df.Point(PLATE_LENGTH,0          ,0),TENSION_LOADS[0])
corner3_load = df.PointSource(V.sub(0), df.Point(0,           PLATE_WIDTH,0),TENSION_LOADS[0])
corner4_load = df.PointSource(V.sub(0), df.Point(PLATE_LENGTH,PLATE_WIDTH,0),TENSION_LOADS[0])

pointsourceb = assemble(Constant(0.0)*v*dx)

corner1_load.apply(pointsourceb)
corner2_load.apply(pointsourceb)
corner3_load.apply(pointsourceb)
corner4_load.apply(pointsourceb)

## Creating the system stiffness K and mass matrix M

k_form = df.inner(stress(du),strain(u_))*df.dx

l_form = df.Constant(1.)*u_[0]*df.dx + pointsourceb
K = df.PETScMatrix()
b = df.PETScVector()

df.assemble_system(k_form, l_form, bc, A_tensor=K, b_tensor=b)

m_form = PLATE_DENSITY*df.dot(du,u_)*df.dx
M = df.PETScMatrix()

assembly = df.assemble(m_form, tensor=M)


# Solving for Modeshapes

eigensolver = df.SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.0
eigensolver.solve(MODE_COUNT)


falpha = lambda x: cos(x)*cosh(x)+1
alpha = lambda n: root(falpha, (2*n+1)*pi/2.)['x'][0]
normal_M =  []
normal_K =  []
lamba = []
Modes = []
Modes_vectors = []
file_results = df.XDMFFile(f"{MAIN_DIR}/modes.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True


def int_zfill(index, total_len): 
    from math import floor, log10
    return str(index).zfill(1+floor(log10(total_len)))

print("Computing {} first eigenvalues...".format(MODE_COUNT))
print("id  | freq(Hz)   | Ang Freq(rad/s) | Lambda")
for i in range(MODE_COUNT):
    # Extract eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    rx[:] /= max(rx[:])

    # 3D eigenfrequency
    freq_3D = df.sqrt(r)/2/np.pi
    lamba.append(r)

    print(f"{i+1:2d}  | {np.abs(freq_3D.real):.3e}  | {np.abs(2*np.pi*freq_3D.real):.3e}       | {np.abs(r.real):.3e}", )

    # Initialize function and assign eigenvector
    eigenmode = df.Function(V,name=f"Eigenvector {int_zfill(i+1, MODE_COUNT)}")
    eigenmode.vector()[:] = rx
    Modes.append(rx[:].reshape((DIM,int(rx[:].shape[0]/DIM)),order='F').T)
    Modes_vectors.append(eigenmode)
    normal_M.append(np.dot(rx[:],M*rx[:]))
    normal_K.append(np.dot(rx[:],K*rx[:]))


    # Function can be plotted as usual:
    rx_func = df.Function(V)
    rx_func.vector()[:] = rx
    file_results.write(eigenmode, 0.)
    file_results.close()

But in getting the following error.

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_992/51715959.py in <module>
    165 # Applying Tensile Loads
    166 
--> 167 corner1_load = df.PointSource(V.sub(0), df.Point(0,           0          ,0),TENSION_LOADS[0])
    168 corner2_load = df.PointSource(V.sub(0), df.Point(PLATE_LENGTH,0          ,0),TENSION_LOADS[0])
    169 corner3_load = df.PointSource(V.sub(0), df.Point(0,           PLATE_WIDTH,0),TENSION_LOADS[0])

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to apply point source to vector.
*** Reason:  The point is outside of the domain.
*** Where:   This error was encountered inside PointSource.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  7f46aeb0b296da5bbb1fb0845822a72ab9b09c55
*** -------------------------------------------------------------------------

This topic is the continuation of the following topic

Any help would be really appreciated.

2

Please note that your code can be massively simplified to illustrate your issue:

# importing libraries for pre-processing

import dolfin as df
import numpy as np
np.set_printoptions(precision=5)

PLATE_LENGTH = 1
WATER_DEPTH = 1

# Plate Data
PLATE_DEPTH = 0.2 * WATER_DEPTH
PLATE_WIDTH = 16 * PLATE_LENGTH
PLATE_THICKNESS = 1e-3 * PLATE_LENGTH

# Meshing Data
#PANEL_HIGH_RES = np.array([20, 160, 6], dtype=int)
PANEL_HIGH_RES = [1, 1, 1]

CORNER1 = df.Point(0, 0, -PLATE_THICKNESS-PLATE_DEPTH)

# corner2 is diagonally opposite to the corner1
CORNER2 = df.Point(PLATE_LENGTH, PLATE_WIDTH, -PLATE_DEPTH)

mesh = df.BoxMesh(CORNER1, CORNER2,
                  PANEL_HIGH_RES[0], PANEL_HIGH_RES[1], PANEL_HIGH_RES[2])

print(mesh.coordinates())
V = df.VectorFunctionSpace(mesh, 'Lagrange', degree=1)
WATER_DENSITY = 1000
GRAVITY = 9.81
TENSION_LOADS = np.array([0.05, 0.15, 1e4]) * \
    (WATER_DENSITY * GRAVITY * WATER_DEPTH**4)
corner1_load = df.PointSource(V.sub(0), df.Point(
    0,           0, 0), TENSION_LOADS[0])
corner2_load = df.PointSource(V.sub(0), df.Point(
    PLATE_LENGTH, 0, 0), TENSION_LOADS[0])
corner3_load = df.PointSource(V.sub(0), df.Point(
    0,           PLATE_WIDTH, 0), TENSION_LOADS[0])
corner4_load = df.PointSource(V.sub(0), df.Point(
    PLATE_LENGTH, PLATE_WIDTH, 0), TENSION_LOADS[0])

If you consider the mesh coordinates, printed by this code, you observe that they are:

[[ 0.     0.    -0.201]
 [ 1.     0.    -0.201]
 [ 0.    16.    -0.201]
 [ 1.    16.    -0.201]
 [ 0.     0.    -0.2  ]
 [ 1.     0.    -0.2  ]
 [ 0.    16.    -0.2  ]
 [ 1.    16.    -0.2  ]]

i.e.
a point with 0 in the z-coordinate is not inside the mesh. This is because of your definition

which you can print:

print(CORNER1.array(), CORNER2.array())

which returns

[ 0.     0.    -0.201] [ 1.  16.  -0.2]
3

Thanks a lot for the reply. Now it is working.

Also, can you tell me how to apply spring boundary conditions for the same corner points of spring stiffness 1e6?

Thanks

4

As I am not familiar with the term spring conditions, it would help if you could write the mathematical formulation of such a condition.

5

Hi,

the presence of springs changes the system stiffness matrix, it cannot be handled by just applying external forces (since the force depends on the displacement of the body) but must be taken into account directly in the variational formulation.

I would suggest replacing your corner springs with distributed springs over a small area as I don’t think you can add “Dirac” terms to a bilinear form (the PointSource functionality but it seems to me that it handles only right-hand sides and not matrices).
So define small subdomains in the corners, mark them with a physical tag (say “1”) and define the ds measure with the corresponding subdomain data.
You can then add to your bilinear form k_form the spring contribution (assuming only vertical springs):

k_form += k_spring*du[2]*u_[2]*ds(1)

where k_spring is a spring stiffness density per unit area
You should not change the rhs.

1 Like