findiff icon indicating copy to clipboard operation
findiff copied to clipboard

Two boundary conditions on the same boundary

Open PetterTh opened this issue 4 years ago • 3 comments

Hi

How would you go about to set two boundary conditions on the same border? For instance in you example with the 1D forced harmonic oscillator with friction to have both u(0)=0 and du/dx (0) = 0.

Petter

PetterTh avatar Jan 22 '20 12:01 PetterTh

Hi, what you describe would be an initial condition problem. Currently findiff only handles boundary value problems. For a given boundary point, you can only specify one condition, otherwise the system would be overdetermined.

maroba avatar Jan 22 '20 17:01 maroba

Ok, thank you for the clarification. But it's not only an issue for initial value problems. Solving the equation for the Euler Beam (https://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_theory#Static_beam_equation) require 4 boundary conditions in total with only two boundaries. Do you think it is something which will be possible to do in the future with this package?

PetterTh avatar Jan 22 '20 20:01 PetterTh

Sorry for the late response. I haven't had time for this project lately. Regarding your question, I must say that the PDE functionality is actually a side feature of this package which actually is only about numerical derivatives. If I have some time, I may go forward into this direction, but please don't count on this in the near future.

maroba avatar Jul 10 '20 12:07 maroba

I believe for cases like solving the beam equation, you create ghost nodes that exist outside of the physical beam and apply the extra boundary conditions to those.

#### SIMPLY SUPPORTED BEAM ####
## Set up the grid
l = 10.0  ## beam length
dx = 0.1
n = int(l//dx + 2)
shape = (n, )
x = np.linspace(-dx, l+dx, shape[0])

## define variables
d = 0.75  ## circular section diameter
q = 10  ## Uniformly distributed load
EpIp = 32e6*np.pi*d**4/64.0 ## Flexural stiffness for 750mm circular element

## Create the PDE
L = EpIp*FinDiff(0, dx, 4)
f = -q*np.ones(shape)

## Apply the boundary conditions for a simply supported beam
bc = BoundaryConditions(shape)
bc[0] = (FinDiff(0,dx,2), 0) ## Enforce zero curvature at left-most ghost point
bc[1] = 0   ## Enforce zero displacement at left support

bc[-1] = (FinDiff(0,dx,2), 0)   ## As above, but for right-most ghost point
bc[-2] = 0  ## As above but for right support

## Solve
pde = PDE(L, f, bc)
y = pde.solve()

## Plot
plt.plot(x,y)
plt.gca().set_xlim(left=0,right=10)
plt.grid()

dangeo84 avatar Sep 18 '23 07:09 dangeo84